diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b5d093e0..bac386104 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ The public API of this library consists of the functions declared in file ## [Unreleased] ### Changed +- Replace internal algorithm for `polygonToCells` with a new version that is more memory-efficient (#785) - Reorganize tests into public / internal. (#762) - Performance enhancement for aarch64, should not affect other platforms (#790, #792) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9ae21f4a1..11198211d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -125,6 +125,7 @@ set(LIB_SOURCE_FILES src/h3lib/include/bbox.h src/h3lib/include/polygon.h src/h3lib/include/polygonAlgos.h + src/h3lib/include/polyfill.h src/h3lib/include/h3Index.h src/h3lib/include/directedEdge.h src/h3lib/include/latLng.h @@ -146,6 +147,7 @@ set(LIB_SOURCE_FILES src/h3lib/lib/coordijk.c src/h3lib/lib/bbox.c src/h3lib/lib/polygon.c + src/h3lib/lib/polyfill.c src/h3lib/lib/h3Index.c src/h3lib/lib/vec2d.c src/h3lib/lib/vec3d.c @@ -192,7 +194,9 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testVertexGraphInternal.c src/apps/testapps/testCompactCells.c src/apps/testapps/testPolygonToCells.c + src/apps/testapps/testPolygonToCellsExperimental.c src/apps/testapps/testPolygonToCellsReported.c + src/apps/testapps/testPolygonToCellsReportedExperimental.c src/apps/testapps/testPentagonIndexes.c src/apps/testapps/testGridDisk.c src/apps/testapps/testGridDiskInternal.c @@ -209,6 +213,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testCellToLatLng.c src/apps/testapps/testCellToCenterChild.c src/apps/testapps/testCellToChildren.c + src/apps/testapps/testCellToBBoxExhaustive.c src/apps/testapps/testCellToChildPos.c src/apps/testapps/testGetIcosahedronFaces.c src/apps/testapps/testLatLng.c @@ -220,6 +225,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testVertexInternal.c src/apps/testapps/testVertexExhaustive.c src/apps/testapps/testPolygonInternal.c + src/apps/testapps/testPolyfillInternal.c src/apps/testapps/testVec2dInternal.c src/apps/testapps/testVec3dInternal.c src/apps/testapps/testDirectedEdge.c @@ -270,6 +276,7 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerInternalAlgos.c src/apps/fuzzers/fuzzerInternalCoordIjk.c src/apps/benchmarks/benchmarkPolygonToCells.c + src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c src/apps/benchmarks/benchmarkPolygon.c src/apps/benchmarks/benchmarkCellsToLinkedMultiPolygon.c src/apps/benchmarks/benchmarkCellToChildren.c @@ -566,6 +573,7 @@ if(BUILD_BENCHMARKS) add_h3_benchmark(benchmarkCellsToLinkedMultiPolygon src/apps/benchmarks/benchmarkCellsToLinkedMultiPolygon.c) add_h3_benchmark(benchmarkCellToChildren src/apps/benchmarks/benchmarkCellToChildren.c) add_h3_benchmark(benchmarkPolygonToCells src/apps/benchmarks/benchmarkPolygonToCells.c) + add_h3_benchmark(benchmarkPolygonToCellsExperimental src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_benchmark(benchmarkPolygon src/apps/benchmarks/benchmarkPolygon.c) endif() diff --git a/CMakeTests.cmake b/CMakeTests.cmake index c8d43a964..a5d57e279 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -186,7 +186,9 @@ add_h3_test(testCellsToLinkedMultiPolygon src/apps/testapps/testCellsToLinkedMul add_h3_test(testH3SetToVertexGraphInternal src/apps/testapps/testH3SetToVertexGraphInternal.c) add_h3_test(testLinkedGeoInternal src/apps/testapps/testLinkedGeoInternal.c) add_h3_test(testPolygonToCells src/apps/testapps/testPolygonToCells.c) +add_h3_test(testPolygonToCellsExperimental src/apps/testapps/testPolygonToCellsExperimental.c) add_h3_test(testPolygonToCellsReported src/apps/testapps/testPolygonToCellsReported.c) +add_h3_test(testPolygonToCellsReportedExperimental src/apps/testapps/testPolygonToCellsReportedExperimental.c) add_h3_test(testVertexGraphInternal src/apps/testapps/testVertexGraphInternal.c) add_h3_test(testDirectedEdge src/apps/testapps/testDirectedEdge.c) add_h3_test(testLatLng src/apps/testapps/testLatLng.c) @@ -195,6 +197,7 @@ add_h3_test(testBBoxInternal src/apps/testapps/testBBoxInternal.c) add_h3_test(testVertex src/apps/testapps/testVertex.c) add_h3_test(testVertexInternal src/apps/testapps/testVertexInternal.c) add_h3_test(testPolygonInternal src/apps/testapps/testPolygonInternal.c) +add_h3_test(testPolyfillInternal src/apps/testapps/testPolyfillInternal.c) add_h3_test(testVec2dInternal src/apps/testapps/testVec2dInternal.c) add_h3_test(testVec3dInternal src/apps/testapps/testVec3dInternal.c) add_h3_test(testCellToLocalIj src/apps/testapps/testCellToLocalIj.c) @@ -223,6 +226,7 @@ add_h3_test(testCellToLocalIjExhaustive src/apps/testapps/testCellToLocalIjExhau add_h3_test(testGridPathCellsExhaustive src/apps/testapps/testGridPathCellsExhaustive.c) add_h3_test(testGridDistanceExhaustive src/apps/testapps/testGridDistanceExhaustive.c) add_h3_test(testH3CellAreaExhaustive src/apps/testapps/testH3CellAreaExhaustive.c) +add_h3_test(testCellToBBoxExhaustive src/apps/testapps/testCellToBBoxExhaustive.c) add_h3_cli_test(testCliCellToLatLng "cellToLatLng -c 8928342e20fffff" "37.5012466151, -122.5003039349") add_h3_cli_test(testCliLatLngToCell "latLngToCell --lat 20 --lng 123 -r 2" "824b9ffffffffff") diff --git a/README.md b/README.md index 8747c0063..a09ff3558 100644 --- a/README.md +++ b/README.md @@ -112,13 +112,13 @@ You can run a faster test suite that excludes the most expensive tests with `mak #### Coverage -You can generate a code coverage report if `lcov` is installed, and if the project was built with the `CMAKE_BUILD_TYPE=Debug` option. +You can generate a code coverage report if `lcov` is installed, and if the project was built with the `CMAKE_BUILD_TYPE=Debug` and `ENABLE_COVERAGE=ON` options. For example, from a clean repository, you could run: ``` mkdir build cd build -cmake -DCMAKE_BUILD_TYPE=Debug .. +cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON .. make make coverage ``` diff --git a/src/apps/applib/include/utility.h b/src/apps/applib/include/utility.h index 0ae5d576e..dcf89e7d1 100644 --- a/src/apps/applib/include/utility.h +++ b/src/apps/applib/include/utility.h @@ -22,6 +22,7 @@ #include +#include "bbox.h" #include "coordijk.h" #include "h3api.h" @@ -46,6 +47,8 @@ void geoPrintNoFmt(const LatLng *p); void geoPrintlnNoFmt(const LatLng *p); void cellBoundaryPrint(const CellBoundary *b); void cellBoundaryPrintln(const CellBoundary *b); +void bboxPrint(const BBox *bbox); +void bboxPrintln(const BBox *bbox); void randomGeo(LatLng *p); diff --git a/src/apps/applib/lib/utility.c b/src/apps/applib/lib/utility.c index 312b1e9ef..15a80cd17 100644 --- a/src/apps/applib/lib/utility.c +++ b/src/apps/applib/lib/utility.c @@ -118,6 +118,18 @@ void cellBoundaryPrintln(const CellBoundary *b) { printf("}\n"); } +void bboxPrint(const BBox *bbox) { + printf( + "bbox {%.9lf, %.9lf, %.9lf, %.9lf}", H3_EXPORT(radsToDegs)(bbox->north), + H3_EXPORT(radsToDegs)(bbox->south), H3_EXPORT(radsToDegs)(bbox->east), + H3_EXPORT(radsToDegs)(bbox->west)); +} + +void bboxPrintln(const BBox *bbox) { + bboxPrint(bbox); + printf("\n"); +} + /** * Apply callback for every unidirectional edge at the given resolution. */ diff --git a/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c b/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c new file mode 100644 index 000000000..b0837c7fc --- /dev/null +++ b/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c @@ -0,0 +1,146 @@ +/* + * Copyright 2017, 2020-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "algos.h" +#include "benchmark.h" +#include "h3api.h" +#include "polyfill.h" + +// Fixtures +LatLng 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}}; +GeoLoop sfGeoLoop; +GeoPolygon sfGeoPolygon; + +LatLng alamedaVerts[] = {{0.6597959342671712, -2.133241848488897}, + {0.6597959348850178, -2.133241848495878}, + {0.6598352639563587, -2.1331688423977755}, + {0.6601346536539207, -2.13270417124178}, + {0.6601594763880223, -2.1326680320633344}, + {0.6601512007732382, -2.1326594176574534}, + {0.6598535076212304, -2.1323049630593562}, + {0.6596565748646488, -2.132069889917591}, + {0.6594645035394391, -2.131843148468039}, + {0.6593438094209757, -2.1316994860539844}, + {0.6591174422311021, -2.131429776816562}, + {0.658849344286881, -2.1311111485483867}, + {0.6588348862079956, -2.1310988536794455}, + {0.6586273138317915, -2.131668420800747}, + {0.6583729538174264, -2.132370426573979}, + {0.6582479206289285, -2.132718691911663}, + {0.6582322393220743, -2.1327614200082317}, + {0.6583003647098981, -2.132837478687196}, + {0.6584457274847966, -2.132827956758973}, + {0.6585526679060995, -2.1330231566043203}, + {0.6587379099516777, -2.1331602726234538}, + {0.6587273684736642, -2.1332676321559063}, + {0.6584638025857692, -2.133305719954319}, + {0.6583545950288919, -2.1334323622944993}, + {0.6584427148370682, -2.1335885223323947}, + {0.6584715236640714, -2.133649780409862}, + {0.6584715242505019, -2.133649780481421}, + {0.658474662092443, -2.1336459234695804}, + {0.6591666596433436, -2.1348354004882926}, + {0.6591809355063646, -2.1348424115474565}, + {0.6593477498700266, -2.1351460576998926}, + {0.6597155087395117, -2.1351049454274}, + {0.6597337410387994, -2.135113899444683}, + {0.6598277083823935, -2.1351065432309517}, + {0.659837290351688, -2.1350919904836627}, + {0.6598391300107502, -2.1350911731005957}, + {0.6598335712627461, -2.1350732321630828}, + {0.6597162034032434, -2.134664026354221}, + {0.6596785831942451, -2.134651647657116}, + {0.6596627824684727, -2.13458880305965}, + {0.6596785832500957, -2.134530719130462}, + {0.6596093592822273, -2.13428052987356}, + {0.6596116166352313, -2.134221493755564}, + {0.6595973199434513, -2.134146270344056}, + {0.6595536764042369, -2.1340805688066653}, + {0.6594611172376618, -2.133753252031165}, + {0.6594829406269346, -2.1337342082305697}, + {0.6594897134102581, -2.1337104032834757}, + {0.6597920983773051, -2.1332343063312775}, + {0.6597959342671712, -2.133241848488897}}; +GeoLoop alamedaGeoLoop; +GeoPolygon alamedaGeoPolygon; + +LatLng southernVerts[] = {{0.6367481147484843, -2.1290865397798906}, + {0.6367481152301953, -2.129086539469222}, + {0.6367550754426818, -2.128887436716856}, + {0.6367816002113981, -2.1273204058681094}, + {0.6380814125349741, -2.127201274803692}, + {0.6388614350074809, -2.12552061082428}, + {0.6393520289210095, -2.124274316938293}, + {0.639524834205869, -2.122168447308359}, + {0.6405714857447717, -2.122083222593005}, + {0.640769478635285, -2.120979885974894}, + {0.6418936996869471, -2.1147667448862255}, + {0.6419094141707652, -2.1146521242709584}, + {0.6269997808948107, -2.1038647304637257}, + {0.6252080524974937, -2.1195521728170457}, + {0.626379700264057, -2.1203708632511162}, + {0.6282200029232767, -2.1210412050690723}, + {0.6283657301211779, -2.1219496416754393}, + {0.6305651783819565, -2.123628532238016}, + {0.6308259852882764, -2.124225549648211}, + {0.6317049665784865, -2.124887756638367}, + {0.6323403882676475, -2.1266205835454053}, + {0.6334397909415498, -2.1277211741619553}, + {0.6367481147484843, -2.1290865397798906}}; +GeoLoop southernGeoLoop; +GeoPolygon southernGeoPolygon; + +BEGIN_BENCHMARKS(); + +sfGeoLoop.numVerts = 6; +sfGeoLoop.verts = sfVerts; +sfGeoPolygon.geoloop = sfGeoLoop; + +alamedaGeoLoop.numVerts = 50; +alamedaGeoLoop.verts = alamedaVerts; +alamedaGeoPolygon.geoloop = alamedaGeoLoop; + +southernGeoLoop.numVerts = 23; +southernGeoLoop.verts = southernVerts; +southernGeoPolygon.geoloop = southernGeoLoop; + +int64_t numHexagons; +H3Index *hexagons; + +BENCHMARK(polygonToCellsSF, 500, { + H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, &numHexagons); + hexagons = calloc(numHexagons, sizeof(H3Index)); + H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0, hexagons); + free(hexagons); +}); + +BENCHMARK(polygonToCellsAlameda, 500, { + H3_EXPORT(maxPolygonToCellsSize)(&alamedaGeoPolygon, 9, 0, &numHexagons); + hexagons = calloc(numHexagons, sizeof(H3Index)); + H3_EXPORT(polygonToCellsExperimental)(&alamedaGeoPolygon, 9, 0, hexagons); + free(hexagons); +}); + +BENCHMARK(polygonToCellsSouthernExpansion, 10, { + H3_EXPORT(maxPolygonToCellsSize)(&southernGeoPolygon, 9, 0, &numHexagons); + hexagons = calloc(numHexagons, sizeof(H3Index)); + H3_EXPORT(polygonToCellsExperimental)(&southernGeoPolygon, 9, 0, hexagons); + free(hexagons); +}); + +END_BENCHMARKS(); diff --git a/src/apps/testapps/testBBoxInternal.c b/src/apps/testapps/testBBoxInternal.c index 8609bbc06..00fff60fa 100644 --- a/src/apps/testapps/testBBoxInternal.c +++ b/src/apps/testapps/testBBoxInternal.c @@ -22,9 +22,10 @@ #include "latLng.h" #include "polygon.h" #include "test.h" +#include "utility.h" -void assertBBox(const GeoLoop *geoloop, const BBox *expected, - const LatLng *inside, const LatLng *outside) { +void assertBBoxFromGeoLoop(const GeoLoop *geoloop, const BBox *expected, + const LatLng *inside, const LatLng *outside) { BBox result; bboxFromGeoLoop(geoloop, &result); @@ -35,14 +36,24 @@ void assertBBox(const GeoLoop *geoloop, const BBox *expected, "Does not contain expected outside point"); } -SUITE(BBoxInternal) { +void assertBBox(const BBox *bbox, const BBox *expected) { + LatLng actualNE = {.lat = bbox->north, .lng = bbox->east}; + LatLng expectedNE = {.lat = expected->north, .lng = expected->east}; + t_assert(geoAlmostEqual(&actualNE, &expectedNE), "NE corner matches"); + + LatLng actualSW = {.lat = bbox->south, .lng = bbox->west}; + LatLng expectedSW = {.lat = expected->south, .lng = expected->west}; + t_assert(geoAlmostEqual(&actualSW, &expectedSW), "SW corner matches"); +} + +SUITE(BBox) { TEST(posLatPosLng) { LatLng verts[] = {{0.8, 0.3}, {0.7, 0.6}, {1.1, 0.7}, {1.0, 0.2}}; const GeoLoop geoloop = {.numVerts = 4, .verts = verts}; const BBox expected = {1.1, 0.7, 0.7, 0.2}; const LatLng inside = {0.9, 0.4}; const LatLng outside = {0.0, 0.0}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(negLatPosLng) { @@ -51,7 +62,7 @@ SUITE(BBoxInternal) { const BBox expected = {-0.1, -0.4, 0.9, 0.6}; const LatLng inside = {-0.3, 0.8}; const LatLng outside = {0.0, 0.0}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(posLatNegLng) { @@ -60,7 +71,7 @@ SUITE(BBoxInternal) { const BBox expected = {1.1, 0.7, -0.8, -1.4}; const LatLng inside = {0.9, -1.0}; const LatLng outside = {0.0, 0.0}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(negLatNegLng) { @@ -70,7 +81,7 @@ SUITE(BBoxInternal) { const BBox expected = {-0.1, -0.4, -1.1, -1.4}; const LatLng inside = {-0.3, -1.2}; const LatLng outside = {0.0, 0.0}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(aroundZeroZero) { @@ -79,7 +90,7 @@ SUITE(BBoxInternal) { const BBox expected = {0.4, -0.4, 0.4, -0.4}; const LatLng inside = {-0.1, -0.1}; const LatLng outside = {1.0, -1.0}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(transmeridian) { @@ -91,7 +102,7 @@ SUITE(BBoxInternal) { const BBox expected = {0.4, -0.4, -M_PI + 0.1, M_PI - 0.1}; const LatLng insideOnMeridian = {-0.1, M_PI}; const LatLng outside = {1.0, M_PI - 0.5}; - assertBBox(&geoloop, &expected, &insideOnMeridian, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &insideOnMeridian, &outside); const LatLng westInside = {0.1, M_PI - 0.05}; t_assert(bboxContains(&expected, &westInside), @@ -117,7 +128,7 @@ SUITE(BBoxInternal) { const BBox expected = {M_PI_2, M_PI_2 - 0.1, 0.8, 0.1}; const LatLng inside = {M_PI_2 - 0.01, 0.4}; const LatLng outside = {M_PI_2, 0.9}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(edgeOnSouthPole) { @@ -129,7 +140,7 @@ SUITE(BBoxInternal) { const BBox expected = {-M_PI_2 + 0.1, -M_PI_2, 0.8, 0.1}; const LatLng inside = {-M_PI_2 + 0.01, 0.4}; const LatLng outside = {-M_PI_2, 0.9}; - assertBBox(&geoloop, &expected, &inside, &outside); + assertBBoxFromGeoLoop(&geoloop, &expected, &inside, &outside); } TEST(containsEdges) { @@ -160,6 +171,99 @@ SUITE(BBoxInternal) { } } + TEST(bboxOverlapsBBox) { + BBox a = {1.0, 0.0, 1.0, 0.0}; + + BBox b1 = {1.0, 0.0, -1.0, -1.5}; + t_assert(!bboxOverlapsBBox(&a, &b1), "no intersection to the west"); + t_assert(!bboxOverlapsBBox(&b1, &a), + "no intersection to the west, reverse"); + + BBox b2 = {1.0, 0.0, 2.0, 1.5}; + t_assert(!bboxOverlapsBBox(&a, &b2), "no intersection to the east"); + t_assert(!bboxOverlapsBBox(&b2, &a), + "no intersection to the east, reverse"); + + BBox b3 = {-1.0, -1.5, 1.0, 0.0}; + t_assert(!bboxOverlapsBBox(&a, &b3), "no intersection to the south"); + t_assert(!bboxOverlapsBBox(&b3, &a), + "no intersection to the south, reverse"); + + BBox b4 = {2.0, 1.5, 1.0, 0.0}; + t_assert(!bboxOverlapsBBox(&a, &b4), "no intersection to the north"); + t_assert(!bboxOverlapsBBox(&b4, &a), + "no intersection to the north, reverse"); + + BBox b5 = {1.0, 0.0, 0.5, -1.5}; + t_assert(bboxOverlapsBBox(&a, &b5), "intersection to the west"); + + BBox b6 = {1.0, 0.0, 2.0, 0.5}; + t_assert(bboxOverlapsBBox(&a, &b6), "intersection to the east"); + + BBox b7 = {0.5, -1.5, 1.0, 0.0}; + t_assert(bboxOverlapsBBox(&a, &b7), "intersection to the south"); + + BBox b8 = {2.0, 0.5, 1.0, 0.0}; + t_assert(bboxOverlapsBBox(&a, &b8), "intersection to the north"); + + BBox b9 = {1.5, -0.5, 1.5, -0.5}; + t_assert(bboxOverlapsBBox(&a, &b9), "intersection, b contains a"); + + BBox b10 = {0.5, 0.25, 0.5, 0.25}; + t_assert(bboxOverlapsBBox(&a, &b10), "intersection, a contains b"); + + BBox b11 = {1.0, 0.0, 1.0, 0.0}; + t_assert(bboxOverlapsBBox(&a, &b11), "intersection, a equals b"); + } + + TEST(bboxOverlapsBBoxTransmeridian) { + BBox a = {1.0, 0.0, -M_PI + 0.5, M_PI - 0.5}; + + BBox b1 = {1.0, 0.0, M_PI - 0.7, M_PI - 0.9}; + t_assert(!bboxOverlapsBBox(&a, &b1), "no intersection to the west"); + t_assert(!bboxOverlapsBBox(&b1, &a), + "no intersection to the west, reverse"); + + BBox b2 = {1.0, 0.0, -M_PI + 0.9, -M_PI + 0.7}; + t_assert(!bboxOverlapsBBox(&a, &b2), "no intersection to the east"); + t_assert(!bboxOverlapsBBox(&b2, &a), "no intersection to the east"); + + BBox b3 = {1.0, 0.0, M_PI - 0.4, M_PI - 0.9}; + t_assert(bboxOverlapsBBox(&a, &b3), "intersection to the west"); + t_assert(bboxOverlapsBBox(&b3, &a), + "intersection to the west, reverse"); + + BBox b4 = {1.0, 0.0, -M_PI + 0.9, -M_PI + 0.4}; + t_assert(bboxOverlapsBBox(&a, &b4), "intersection to the east"); + t_assert(bboxOverlapsBBox(&b4, &a), + "intersection to the east, reverse"); + + BBox b5 = {1.0, 0.0, -M_PI + 0.4, M_PI - 0.4}; + t_assert(bboxOverlapsBBox(&a, &b5), "intersection, a contains b"); + t_assert(bboxOverlapsBBox(&b5, &a), + "intersection, a contains b, reverse"); + + BBox b6 = {1.0, 0.0, -M_PI + 0.6, M_PI - 0.6}; + t_assert(bboxOverlapsBBox(&a, &b6), "intersection, b contains a"); + t_assert(bboxOverlapsBBox(&b6, &a), + "intersection, b contains a, reverse"); + + BBox b7 = {1.0, 0.0, -M_PI + 0.5, M_PI - 0.5}; + t_assert(bboxOverlapsBBox(&a, &b7), "intersection, a equals b"); + + BBox b8 = {1.0, 0.0, -M_PI + 0.9, M_PI - 0.4}; + t_assert(bboxOverlapsBBox(&a, &b8), + "intersection, transmeridian to the east"); + t_assert(bboxOverlapsBBox(&b8, &a), + "intersection, transmeridian to the east, reverse"); + + BBox b9 = {1.0, 0.0, -M_PI + 0.4, M_PI - 0.9}; + t_assert(bboxOverlapsBBox(&a, &b9), + "intersection, transmeridian to the west"); + t_assert(bboxOverlapsBBox(&b9, &a), + "intersection, transmeridian to the west, reverse"); + } + TEST(bboxCenterBasicQuandrants) { LatLng center; @@ -274,4 +378,60 @@ SUITE(BBoxInternal) { E_RES_DOMAIN, "lineHexEstimate of invalid resolution fails"); } + + TEST(scaleBBox_noop) { + BBox bbox = {1.0, 0.0, 1.0, 0.0}; + BBox expected = {1.0, 0.0, 1.0, 0.0}; + scaleBBox(&bbox, 1); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_basicGrow) { + BBox bbox = {1.0, 0.0, 1.0, 0.0}; + BBox expected = {1.5, -0.5, 1.5, -0.5}; + scaleBBox(&bbox, 2); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_basicShrink) { + BBox bbox = {1.0, 0.0, 1.0, 0.0}; + BBox expected = {0.75, 0.25, 0.75, 0.25}; + scaleBBox(&bbox, 0.5); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_clampNorthSouth) { + BBox bbox = {M_PI_2 * 0.9, -M_PI_2 * 0.9, 1.0, 0.0}; + BBox expected = {M_PI_2, -M_PI_2, 1.5, -0.5}; + scaleBBox(&bbox, 2); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_clampEastPos) { + BBox bbox = {1.0, 0.0, M_PI - 0.1, M_PI - 1.1}; + BBox expected = {1.5, -0.5, -M_PI + 0.4, M_PI - 1.6}; + scaleBBox(&bbox, 2); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_clampEastNeg) { + BBox bbox = {1.5, -0.5, -M_PI + 0.4, M_PI - 1.6}; + BBox expected = {1.0, 0.0, M_PI - 0.1, M_PI - 1.1}; + scaleBBox(&bbox, 0.5); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_clampWestPos) { + BBox bbox = {1.0, 0.0, -M_PI + 0.9, M_PI - 0.1}; + BBox expected = {0.75, 0.25, -M_PI + 0.65, -M_PI + 0.15}; + scaleBBox(&bbox, 0.5); + assertBBox(&bbox, &expected); + } + + TEST(scaleBBox_clampWestNeg) { + BBox bbox = {0.75, 0.25, -M_PI + 0.65, -M_PI + 0.15}; + BBox expected = {1.0, 0.0, -M_PI + 0.9, M_PI - 0.1}; + scaleBBox(&bbox, 2); + assertBBox(&bbox, &expected); + } } diff --git a/src/apps/testapps/testCellToBBoxExhaustive.c b/src/apps/testapps/testCellToBBoxExhaustive.c new file mode 100644 index 000000000..bfbbeb4a3 --- /dev/null +++ b/src/apps/testapps/testCellToBBoxExhaustive.c @@ -0,0 +1,100 @@ +/* + * Copyright 2023 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file + * @brief tests H3 function `cellToBBox` + * + * usage: `testCellToBBoxExhaustive` + */ + +#include +#include +#include + +#include "bbox.h" +#include "h3Index.h" +#include "h3api.h" +#include "polyfill.h" +#include "polygon.h" +#include "test.h" +#include "utility.h" + +static void cellBBox_assertions(H3Index h3) { + BBox bbox; + t_assertSuccess(cellToBBox(h3, &bbox, false)); + + CellBoundary verts; + t_assertSuccess(H3_EXPORT(cellToBoundary)(h3, &verts)); + for (int j = 0; j < verts.numVerts; j++) { + if (!bboxContains(&bbox, &verts.verts[j])) { + printf("cell: "); + h3Println(h3); + bboxPrintln(&bbox); + geoPrintln(&verts.verts[j]); + } + t_assert(bboxContains(&bbox, &verts.verts[j]), + "BBox contains cell vertex"); + } +} + +static void childBBox_assertions(H3Index h3) { + int parentRes = H3_GET_RESOLUTION(h3); + + BBox bbox; + t_assertSuccess(cellToBBox(h3, &bbox, true)); + + for (int resolutionOffset = 0; resolutionOffset < 5; resolutionOffset++) { + // Test whether all verts of all children are inside the bbox + int childRes = parentRes + resolutionOffset; + int64_t numChildren; + t_assertSuccess( + H3_EXPORT(cellToChildrenSize)(h3, childRes, &numChildren)); + + H3Index *children = calloc(numChildren, sizeof(H3Index)); + t_assertSuccess(H3_EXPORT(cellToChildren)(h3, childRes, children)); + + for (int64_t i = 0; i < numChildren; i++) { + CellBoundary childVerts; + t_assertSuccess( + H3_EXPORT(cellToBoundary)(children[i], &childVerts)); + for (int j = 0; j < childVerts.numVerts; j++) { + if (!bboxContains(&bbox, &childVerts.verts[j])) { + printf("parent: "); + h3Println(h3); + bboxPrintln(&bbox); + h3Println(children[i]); + geoPrintln(&childVerts.verts[j]); + } + t_assert(bboxContains(&bbox, &childVerts.verts[j]), + "BBox contains child vertex"); + } + } + + free(children); + } +} + +SUITE(cellToBBox) { + TEST(cellBBox_correctness) { + iterateAllIndexesAtRes(0, cellBBox_assertions); + iterateAllIndexesAtRes(1, cellBBox_assertions); + iterateAllIndexesAtRes(2, cellBBox_assertions); + } + TEST(childBBox_correctness) { + iterateAllIndexesAtRes(0, childBBox_assertions); + iterateAllIndexesAtRes(1, childBBox_assertions); + iterateAllIndexesAtRes(2, childBBox_assertions); + } +} diff --git a/src/apps/testapps/testH3Memory.c b/src/apps/testapps/testH3Memory.c index 3bed50937..96627b543 100644 --- a/src/apps/testapps/testH3Memory.c +++ b/src/apps/testapps/testH3Memory.c @@ -25,6 +25,7 @@ #include "h3Index.h" #include "h3api.h" #include "latLng.h" +#include "polyfill.h" #include "test.h" #include "utility.h" @@ -223,4 +224,35 @@ SUITE(h3Memory) { t_assert(actualNumIndexes == 1253, "got expected polygonToCells size"); free(hexagons); } + + TEST(polygonToCellsExperimental) { + sfGeoPolygon.geoloop = sfGeoLoop; + sfGeoPolygon.numHoles = 0; + + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + resetMemoryCounters(0); + failAlloc = true; + H3Error err = H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0, + hexagons); + t_assert(err == E_MEMORY_ALLOC, + "polygonToCellsExperimental failed (1)"); + t_assert(actualAllocCalls == 1, "alloc called once"); + t_assert(actualFreeCalls == 0, "free not called"); + + resetMemoryCounters(1); + err = H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0, + hexagons); + t_assert(err == E_SUCCESS, "polygonToCellsExperimental succeeded (1)"); + t_assert(actualAllocCalls == 1, "alloc called one time"); + t_assert(actualFreeCalls == 1, "free called one time"); + + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + t_assert(actualNumIndexes == 1253, + "got expected polygonToCellsExperimental size"); + free(hexagons); + } } diff --git a/src/apps/testapps/testPolyfillInternal.c b/src/apps/testapps/testPolyfillInternal.c new file mode 100644 index 000000000..db393cc30 --- /dev/null +++ b/src/apps/testapps/testPolyfillInternal.c @@ -0,0 +1,172 @@ +/* + * Copyright 2017-2018, 2020-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "bbox.h" +#include "h3Index.h" +#include "h3api.h" +#include "latLng.h" +#include "polyfill.h" +#include "test.h" +#include "utility.h" + +// Fixtures +static GeoPolygon sfGeoPolygon = { + .geoloop = {.numVerts = 6, + .verts = (LatLng[]){{0.659966917655, -2.1364398519396}, + {0.6595011102219, -2.1359434279405}, + {0.6583348114025, -2.1354884206045}, + {0.6581220034068, -2.1382437718946}, + {0.6594479998527, -2.1384597563896}, + {0.6599990002976, -2.1376771158464}}}, + .numHoles = 0}; + +SUITE(polyfillInternal) { + TEST(iterInitPolygonCompact_errors) { + IterCellsPolygonCompact iter; + + iter = iterInitPolygonCompact(&sfGeoPolygon, -1, 0); + t_assert(iter.error == E_RES_DOMAIN, + "Got expected error for invalid res"); + t_assert(iter.cell == H3_NULL, "Got null output for invalid res"); + + iter = iterInitPolygonCompact(&sfGeoPolygon, 16, 0); + t_assert(iter.error == E_RES_DOMAIN, + "Got expected error for invalid res"); + t_assert(iter.cell == H3_NULL, "Got null output for invalid res"); + + iter = iterInitPolygonCompact(&sfGeoPolygon, 9, 42); + t_assert(iter.error == E_OPTION_INVALID, + "Got expected error for invalid flags"); + t_assert(iter.cell == H3_NULL, "Got null output for invalid flags"); + } + + TEST(iterStepPolygonCompact_errors) { + IterCellsPolygonCompact iter; + H3Index cell; + + iter = iterInitPolygonCompact(&sfGeoPolygon, 9, 0); + + // Give the iterator a cell with a bad base cell + cell = 0x85283473fffffff; + H3_SET_BASE_CELL(cell, 123); + iter.cell = cell; + + iterStepPolygonCompact(&iter); + t_assert(iter.error == E_CELL_INVALID, + "Got expected error for invalid cell"); + t_assert(iter.cell == H3_NULL, "Got null output for invalid cell"); + + iter = iterInitPolygonCompact(&sfGeoPolygon, 9, 0); + + // Give the iterator a cell that's too fine for a child check, + // and a target resolution that allows this to run. This cell has + // to be inside the polygon to reach the error. + cell = 0x8f283080dcb019a; + iter.cell = cell; + iter._res = 42; + + iterStepPolygonCompact(&iter); + t_assert(iter.error == E_RES_DOMAIN, + "Got expected error for too-fine cell"); + t_assert(iter.cell == H3_NULL, "Got null output for invalid cell"); + } + + TEST(iterDestroyPolygonCompact) { + IterCellsPolygonCompact iter = + iterInitPolygonCompact(&sfGeoPolygon, 9, 0); + + iterDestroyPolygonCompact(&iter); + t_assert(iter.error == E_SUCCESS, "Got success for destroyed iterator"); + t_assert(iter.cell == H3_NULL, + "Got null output for destroyed iterator"); + + for (int i = 0; i < 3; i++) { + iterStepPolygonCompact(&iter); + t_assert(iter.cell == H3_NULL, + "Got null output for destroyed iterator"); + } + } + + TEST(iterDestroyPolygon) { + IterCellsPolygon iter = iterInitPolygon(&sfGeoPolygon, 9, 0); + + iterDestroyPolygon(&iter); + t_assert(iter.error == E_SUCCESS, "Got success for destroyed iterator"); + t_assert(iter.cell == H3_NULL, + "Got null output for destroyed iterator"); + + for (int i = 0; i < 3; i++) { + iterStepPolygon(&iter); + t_assert(iter.cell == H3_NULL, + "Got null output for destroyed iterator"); + } + } + + TEST(cellToBBox_noScale) { + // arbitrary cell + H3Index cell = 0x85283473fffffff; + BBox bbox; + t_assertSuccess(cellToBBox(cell, &bbox, false)); + + double cellArea; + t_assertSuccess(H3_EXPORT(cellAreaRads2)(cell, &cellArea)); + double bboxArea = bboxWidthRads(&bbox) * bboxHeightRads(&bbox); + double ratio = bboxArea / cellArea; + + CellBoundary verts; + t_assertSuccess(H3_EXPORT(cellToBoundary)(cell, &verts)); + + t_assert(ratio < 3 && ratio > 1, + "Got reasonable area ratio between cell and bbox"); + } + + TEST(cellToBBox_boundaryError) { + // arbitrary cell + H3Index cell = 0x85283473fffffff; + H3_SET_BASE_CELL(cell, 123); + + BBox bbox; + t_assert(cellToBBox(cell, &bbox, false) == E_CELL_INVALID, + "Got expected error for cell with invalid base cell"); + } + + TEST(cellToBBox_res0boundaryError) { + // arbitrary res 0 cell + H3Index cell = 0x8001fffffffffff; + H3_SET_BASE_CELL(cell, 123); + + BBox bbox; + t_assert(cellToBBox(cell, &bbox, false) == E_CELL_INVALID, + "Got expected error for res 0 cell with invalid base cell"); + } + + TEST(baseCellNumToCell) { + for (int i = 0; i < NUM_BASE_CELLS; i++) { + H3Index cell = baseCellNumToCell(i); + t_assert(H3_EXPORT(isValidCell)(cell), "Cell is valid"); + t_assert( + H3_GET_BASE_CELL(cell) == i && H3_GET_RESOLUTION(cell) == 0, + "Cell has correct base number and res"); + } + } + + TEST(baseCellNumToCell_boundaryErrors) { + t_assert(baseCellNumToCell(-1) == H3_NULL, + "Expected null for less than 0"); + t_assert(baseCellNumToCell(NUM_BASE_CELLS) == H3_NULL, + "Expected null for positive out of bounds"); + } +} diff --git a/src/apps/testapps/testPolygonInternal.c b/src/apps/testapps/testPolygonInternal.c index 3994c9902..4b0f80286 100644 --- a/src/apps/testapps/testPolygonInternal.c +++ b/src/apps/testapps/testPolygonInternal.c @@ -625,4 +625,233 @@ SUITE(polygonInternal) { H3_EXPORT(destroyLinkedMultiPolygon)(&polygon); } + + TEST(lineIntersectsLine) { + LatLng lines1[4] = {{0, 0}, {1, 1}, {0, 1}, {1, 0}}; + t_assert( + lineIntersectsLine(&lines1[0], &lines1[1], &lines1[2], &lines1[3]), + "diagonal intersection"); + + LatLng lines2[4] = {{1, 1}, {0, 0}, {1, 0}, {0, 1}}; + t_assert( + lineIntersectsLine(&lines2[0], &lines2[1], &lines2[2], &lines2[3]), + "diagonal intersection, reverse vertexes"); + + LatLng lines3[4] = {{0.5, 0}, {0.5, 1}, {0, 0.5}, {1, 0.5}}; + t_assert( + lineIntersectsLine(&lines3[0], &lines3[1], &lines3[2], &lines3[3]), + "horizontal/vertical intersection"); + + LatLng lines4[4] = {{0.5, 1}, {0.5, 0}, {1, 0.5}, {0, 0.5}}; + t_assert( + lineIntersectsLine(&lines4[0], &lines4[1], &lines4[2], &lines4[3]), + "horizontal/vertical intersection, reverse vertexes"); + + LatLng lines5[4] = {{0, 0}, {0.4, 0.4}, {0, 1}, {1, 0}}; + t_assert( + !lineIntersectsLine(&lines5[0], &lines5[1], &lines5[2], &lines5[3]), + "diagonal non-intersection, below"); + + LatLng lines6[4] = {{0.6, 0.6}, {1, 1}, {0, 1}, {1, 0}}; + t_assert( + !lineIntersectsLine(&lines6[0], &lines6[1], &lines6[2], &lines6[3]), + "diagonal non-intersection, above"); + + LatLng lines7[4] = {{0.5, 0}, {0.5, 1}, {0, 0.5}, {0.4, 0.5}}; + t_assert( + !lineIntersectsLine(&lines7[0], &lines7[1], &lines7[2], &lines7[3]), + "horizontal/vertical non-intersection, below"); + + LatLng lines8[4] = {{0.5, 0}, {0.5, 1}, {0.6, 0.5}, {1, 0.5}}; + t_assert( + !lineIntersectsLine(&lines8[0], &lines8[1], &lines8[2], &lines8[3]), + "horizontal/vertical non-intersection, above"); + + LatLng lines9[4] = {{0.5, 0}, {0.5, 0.4}, {0, 0.5}, {1, 0.5}}; + t_assert( + !lineIntersectsLine(&lines9[0], &lines9[1], &lines9[2], &lines9[3]), + "horizontal/vertical non-intersection, left"); + + LatLng lines10[4] = {{0.5, 0.6}, {0.5, 1}, {0, 0.5}, {1, 0.5}}; + t_assert(!lineIntersectsLine(&lines10[0], &lines10[1], &lines10[2], + &lines10[3]), + "horizontal/vertical non-intersection, right"); + } + + TEST(cellBoundaryInsidePolygon_inside) { + LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; + 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.6, 0.6}, {0.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}}}; + BBox boundaryBBox = {0.6, 0.4, 0.6, 0.4}; + + t_assert(cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "simple containment is inside"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_insideTransmeridianWest) { + LatLng verts[] = {{0, M_PI - 0.5}, + {0, -M_PI + 0.5}, + {1, -M_PI + 0.5}, + {1, M_PI - 0.5}}; + 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.6, M_PI - 0.1}, + {0.6, M_PI - 0.2}, + {0.4, M_PI - 0.2}, + {0.4, M_PI - 0.1}}}; + BBox boundaryBBox = {0.6, 0.4, 0.6, 0.4}; + + t_assert(cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "simple containment is inside, west side of transmeridian"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_insideTransmeridianEast) { + LatLng verts[] = {{0, M_PI - 0.5}, + {0, -M_PI + 0.5}, + {1, -M_PI + 0.5}, + {1, M_PI - 0.5}}; + 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.6, -M_PI + 0.4}, + {0.6, -M_PI + 0.2}, + {0.4, -M_PI + 0.2}, + {0.4, -M_PI + 0.4}}}; + BBox boundaryBBox = {0.6, 0.4, 0.6, 0.4}; + + t_assert(cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "simple containment is inside, east side of transmeridian"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_insideWithHole) { + LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; + GeoLoop geoloop = {.numVerts = 4, .verts = verts}; + + LatLng holeVerts[] = {{0.3, 0.3}, {0.3, 0.1}, {0.1, 0.1}, {0.1, 0.3}}; + 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.6, 0.6}, {0.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}}}; + BBox boundaryBBox = {0.6, 0.4, 0.6, 0.4}; + + t_assert(cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "simple containment is inside, with hole"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_notInside) { + LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; + 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 = {{1.6, 1.6}, {1.6, 1.4}, {1.4, 1.4}, {1.4, 1.6}}}; + BBox boundaryBBox = {1.6, 1.4, 1.6, 1.4}; + + t_assert(!cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "fully outside is not inside"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_notInsideIntersect) { + LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; + 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.6, 0.6}, {1.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}}}; + BBox boundaryBBox = {1.6, 0.4, 0.6, 0.4}; + + t_assert(!cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "intersecting polygon is not inside"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_notInsideIntersectHole) { + LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; + GeoLoop geoloop = {.numVerts = 4, .verts = verts}; + + LatLng holeVerts[] = {{0.3, 0.3}, {0.5, 0.5}, {0.1, 0.1}, {0.1, 0.3}}; + 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.6, 0.6}, {0.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}}}; + BBox boundaryBBox = {0.6, 0.4, 0.6, 0.4}; + + t_assert(!cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "not inside with hole intersection"); + free(bboxes); + } + + TEST(cellBoundaryInsidePolygon_notInsideWithinHole) { + LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}; + GeoLoop geoloop = {.numVerts = 4, .verts = verts}; + + LatLng holeVerts[] = {{0.9, 0.9}, {0.9, 0.1}, {0.1, 0.1}, {0.1, 0.9}}; + 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.6, 0.6}, {0.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}}}; + BBox boundaryBBox = {0.6, 0.4, 0.6, 0.4}; + + t_assert(!cellBoundaryInsidePolygon(&polygon, bboxes, &boundary, + &boundaryBBox), + "not inside when within hole"); + free(bboxes); + } } diff --git a/src/apps/testapps/testPolygonToCells.c b/src/apps/testapps/testPolygonToCells.c index 65050463c..22a0cdecf 100644 --- a/src/apps/testapps/testPolygonToCells.c +++ b/src/apps/testapps/testPolygonToCells.c @@ -158,6 +158,10 @@ SUITE(polygonToCells) { lineGeoPolygon.geoloop = lineGeoLoop; lineGeoPolygon.numHoles = 0; + // -------------------------------------------- + // maxPolygonToCellsSize + // -------------------------------------------- + TEST(maxPolygonToCellsSize) { int64_t numHexagons; t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, @@ -175,6 +179,35 @@ SUITE(polygonToCells) { "got expected max polygonToCells size (empty)"); } + TEST(maxPolygonToCellsSizeInvalid) { + int64_t numHexagons; + t_assert( + H3_EXPORT(maxPolygonToCellsSize)(&invalidGeoPolygon, 9, 0, + &numHexagons) == E_FAILED, + "Cannot determine cell size to invalid geo polygon with Infinity"); + t_assert(H3_EXPORT(maxPolygonToCellsSize)(&invalid2GeoPolygon, 9, 0, + &numHexagons) == E_FAILED, + "Cannot determine cell size to invalid geo polygon with NaNs"); + } + + TEST(maxPolygonToCellsSizePoint) { + int64_t numHexagons; + t_assert(H3_EXPORT(maxPolygonToCellsSize)(&pointGeoPolygon, 9, 0, + &numHexagons) == E_FAILED, + "Cannot estimate for single point"); + } + + TEST(maxPolygonToCellsSizeLine) { + int64_t numHexagons; + t_assert(H3_EXPORT(maxPolygonToCellsSize)(&lineGeoPolygon, 9, 0, + &numHexagons) == E_FAILED, + "Cannot estimate for straight line"); + } + + // -------------------------------------------- + // polygonToCells + // -------------------------------------------- + TEST(polygonToCells) { int64_t numHexagons; t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, @@ -445,6 +478,18 @@ SUITE(polygonToCells) { free(hexagons); } + TEST(polygonToCellsInvalidPolygon) { + // Chosen arbitrarily, polygonToCells should error out before this is an + // issue. + int64_t numHexagons = 0; + + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + t_assert(H3_EXPORT(polygonToCells)(&invalidGeoPolygon, 9, 0, + hexagons) == E_FAILED, + "Invalid geo polygon cannot be evaluated"); + free(hexagons); + } + TEST(fillIndex) { iterateAllIndexesAtRes(0, fillIndex_assertions); iterateAllIndexesAtRes(1, fillIndex_assertions); @@ -468,39 +513,4 @@ SUITE(polygonToCells) { free(found); free(search); } - - TEST(polygonToCellsInvalid) { - int64_t numHexagons; - t_assert( - H3_EXPORT(maxPolygonToCellsSize)(&invalidGeoPolygon, 9, 0, - &numHexagons) == E_FAILED, - "Cannot determine cell size to invalid geo polygon with Infinity"); - t_assert(H3_EXPORT(maxPolygonToCellsSize)(&invalid2GeoPolygon, 9, 0, - &numHexagons) == E_FAILED, - "Cannot determine cell size to invalid geo polygon with NaNs"); - - // Chosen arbitrarily, polygonToCells should error out before this is an - // issue. - numHexagons = 0; - - H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); - t_assert(H3_EXPORT(polygonToCells)(&invalidGeoPolygon, 9, 0, - hexagons) == E_FAILED, - "Invalid geo polygon cannot be evaluated"); - free(hexagons); - } - - TEST(polygonToCellsPoint) { - int64_t numHexagons; - t_assert(H3_EXPORT(maxPolygonToCellsSize)(&pointGeoPolygon, 9, 0, - &numHexagons) == E_FAILED, - "Cannot estimate for single point"); - } - - TEST(polygonToCellsLine) { - int64_t numHexagons; - t_assert(H3_EXPORT(maxPolygonToCellsSize)(&lineGeoPolygon, 9, 0, - &numHexagons) == E_FAILED, - "Cannot estimate for straight line"); - } } diff --git a/src/apps/testapps/testPolygonToCellsExperimental.c b/src/apps/testapps/testPolygonToCellsExperimental.c new file mode 100644 index 000000000..f389456ab --- /dev/null +++ b/src/apps/testapps/testPolygonToCellsExperimental.c @@ -0,0 +1,438 @@ +/* + * Copyright 2017-2018, 2020-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include + +#include "constants.h" +#include "h3Index.h" +#include "latLng.h" +#include "polyfill.h" +#include "test.h" +#include "utility.h" + +// Fixtures +static LatLng 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}}; +static GeoLoop sfGeoLoop = {.numVerts = 6, .verts = sfVerts}; +static GeoPolygon sfGeoPolygon; + +static LatLng holeVerts[] = {{0.6595072188743, -2.1371053983433}, + {0.6591482046471, -2.1373141048153}, + {0.6592295020837, -2.1365222838402}}; +static GeoLoop holeGeoLoop = {.numVerts = 3, .verts = holeVerts}; +static GeoPolygon holeGeoPolygon; + +static LatLng emptyVerts[] = {{0.659966917655, -2.1364398519394}, + {0.659966917656, -2.1364398519395}, + {0.659966917657, -2.1364398519396}}; +static GeoLoop emptyGeoLoop = {.numVerts = 3, .verts = emptyVerts}; +static GeoPolygon emptyGeoPolygon; + +static LatLng invalidVerts[] = {{INFINITY, INFINITY}, {-INFINITY, -INFINITY}}; +static GeoLoop invalidGeoLoop = {.numVerts = 2, .verts = invalidVerts}; +static GeoPolygon invalidGeoPolygon; + +static LatLng invalid2Verts[] = {{NAN, NAN}, {-NAN, -NAN}}; +static GeoLoop invalid2GeoLoop = {.numVerts = 2, .verts = invalid2Verts}; +static GeoPolygon invalid2GeoPolygon; + +static LatLng pointVerts[] = {{0, 0}}; +static GeoLoop pointGeoLoop = {.numVerts = 1, .verts = pointVerts}; +static GeoPolygon pointGeoPolygon; + +static LatLng lineVerts[] = {{0, 0}, {1, 0}}; +static GeoLoop lineGeoLoop = {.numVerts = 2, .verts = lineVerts}; +static GeoPolygon lineGeoPolygon; + +/** + * Return true if the cell crosses the meridian. + */ +static bool isTransmeridianCell(H3Index h) { + CellBoundary bndry; + H3_EXPORT(cellToBoundary)(h, &bndry); + + double minLng = M_PI, maxLng = -M_PI; + for (int i = 0; i < bndry.numVerts; i++) { + if (bndry.verts[i].lng < minLng) minLng = bndry.verts[i].lng; + if (bndry.verts[i].lng > maxLng) maxLng = bndry.verts[i].lng; + } + + return maxLng - minLng > M_PI - (M_PI / 4); +} + +static void fillIndex_assertions(H3Index h) { + if (isTransmeridianCell(h)) { + // TODO: these do not work correctly + return; + } + + int currentRes = H3_EXPORT(getResolution)(h); + // TODO: Not testing more than one depth because the assertions fail. + for (int nextRes = currentRes; nextRes <= currentRes + 1; nextRes++) { + CellBoundary bndry; + H3_EXPORT(cellToBoundary)(h, &bndry); + GeoPolygon polygon = { + .geoloop = {.numVerts = bndry.numVerts, .verts = bndry.verts}, + .numHoles = 0, + .holes = 0}; + + int64_t polygonToCellsSize; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&polygon, nextRes, 0, + &polygonToCellsSize)); + H3Index *polygonToCellsOut = + calloc(polygonToCellsSize, sizeof(H3Index)); + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &polygon, nextRes, 0, polygonToCellsOut)); + + int64_t polygonToCellsCount = + countNonNullIndexes(polygonToCellsOut, polygonToCellsSize); + + int64_t childrenSize; + t_assertSuccess( + H3_EXPORT(cellToChildrenSize)(h, nextRes, &childrenSize)); + H3Index *children = calloc(childrenSize, sizeof(H3Index)); + t_assertSuccess(H3_EXPORT(cellToChildren)(h, nextRes, children)); + + int64_t cellToChildrenCount = + countNonNullIndexes(children, childrenSize); + + t_assert(polygonToCellsCount == cellToChildrenCount, + "PolygonToCells count matches cellToChildren count"); + + for (int i = 0; i < childrenSize; i++) { + bool found = false; + if (children[i] == H3_NULL) continue; + for (int j = 0; j < polygonToCellsSize; j++) { + if (polygonToCellsOut[j] == children[i]) { + found = true; + break; + } + } + t_assert( + found, + "All indexes match between polygonToCells and cellToChildren"); + } + + free(polygonToCellsOut); + free(children); + } +} + +SUITE(polygonToCells) { + sfGeoPolygon.geoloop = sfGeoLoop; + sfGeoPolygon.numHoles = 0; + + holeGeoPolygon.geoloop = sfGeoLoop; + holeGeoPolygon.numHoles = 1; + holeGeoPolygon.holes = &holeGeoLoop; + + emptyGeoPolygon.geoloop = emptyGeoLoop; + emptyGeoPolygon.numHoles = 0; + + invalidGeoPolygon.geoloop = invalidGeoLoop; + invalidGeoPolygon.numHoles = 0; + + invalid2GeoPolygon.geoloop = invalid2GeoLoop; + invalid2GeoPolygon.numHoles = 0; + + pointGeoPolygon.geoloop = pointGeoLoop; + pointGeoPolygon.numHoles = 0; + + lineGeoPolygon.geoloop = lineGeoLoop; + lineGeoPolygon.numHoles = 0; + + TEST(polygonToCells) { + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, + 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 1253, "got expected polygonToCells size"); + free(hexagons); + } + + TEST(polygonToCellsHole) { + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&holeGeoPolygon, 9, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&holeGeoPolygon, + 9, 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 1214, + "got expected polygonToCells size (hole)"); + free(hexagons); + } + + TEST(polygonToCellsEmpty) { + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&emptyGeoPolygon, 9, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&emptyGeoPolygon, + 9, 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 0, + "got expected polygonToCells size (empty)"); + free(hexagons); + } + + TEST(polygonToCellsExact) { + LatLng somewhere = {1, 2}; + H3Index origin; + t_assertSuccess(H3_EXPORT(latLngToCell)(&somewhere, 9, &origin)); + CellBoundary boundary; + H3_EXPORT(cellToBoundary)(origin, &boundary); + + LatLng *verts = calloc(boundary.numVerts + 1, sizeof(LatLng)); + for (int i = 0; i < boundary.numVerts; i++) { + verts[i] = boundary.verts[i]; + } + verts[boundary.numVerts] = boundary.verts[0]; + + GeoLoop someGeoLoop; + someGeoLoop.numVerts = boundary.numVerts + 1; + someGeoLoop.verts = verts; + GeoPolygon someHexagon; + someHexagon.geoloop = someGeoLoop; + someHexagon.numHoles = 0; + + int64_t numHexagons; + t_assertSuccess( + H3_EXPORT(maxPolygonToCellsSize)(&someHexagon, 9, 0, &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&someHexagon, 9, + 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 1, "got expected polygonToCells size (1)"); + free(hexagons); + free(verts); + } + + TEST(polygonToCellsTransmeridian) { + LatLng primeMeridianVerts[] = { + {0.01, 0.01}, {0.01, -0.01}, {-0.01, -0.01}, {-0.01, 0.01}}; + GeoLoop primeMeridianGeoLoop = {.numVerts = 4, + .verts = primeMeridianVerts}; + GeoPolygon primeMeridianGeoPolygon = {.geoloop = primeMeridianGeoLoop, + .numHoles = 0}; + + LatLng 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}}; + GeoLoop transMeridianGeoLoop = {.numVerts = 4, + .verts = transMeridianVerts}; + GeoPolygon transMeridianGeoPolygon = {.geoloop = transMeridianGeoLoop, + .numHoles = 0}; + + LatLng 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}}; + GeoLoop transMeridianHoleGeoLoop = {.numVerts = 4, + .verts = transMeridianHoleVerts}; + GeoPolygon transMeridianHoleGeoPolygon = { + .geoloop = transMeridianGeoLoop, + .numHoles = 1, + .holes = &transMeridianHoleGeoLoop}; + GeoPolygon transMeridianFilledHoleGeoPolygon = { + .geoloop = transMeridianHoleGeoLoop, .numHoles = 0}; + + int64_t expectedSize; + + // Prime meridian case + expectedSize = 4228; + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)( + &primeMeridianGeoPolygon, 7, 0, &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &primeMeridianGeoPolygon, 7, 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == expectedSize, + "got expected polygonToCells size (prime meridian)"); + + // Transmeridian case + // This doesn't exactly match the prime meridian count because of slight + // differences in hex size and grid offset between the two cases + expectedSize = 4238; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)( + &transMeridianGeoPolygon, 7, 0, &numHexagons)); + H3Index *hexagonsTM = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &transMeridianGeoPolygon, 7, 0, hexagonsTM)); + actualNumIndexes = countNonNullIndexes(hexagonsTM, numHexagons); + + t_assert(actualNumIndexes == expectedSize, + "got expected polygonToCells size (transmeridian)"); + + // Transmeridian filled hole case -- only needed for calculating hole + // size + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)( + &transMeridianFilledHoleGeoPolygon, 7, 0, &numHexagons)); + H3Index *hexagonsTMFH = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &transMeridianFilledHoleGeoPolygon, 7, 0, hexagonsTMFH)); + int64_t actualNumHoleIndexes = + countNonNullIndexes(hexagonsTMFH, numHexagons); + + // Transmeridian hole case + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)( + &transMeridianHoleGeoPolygon, 7, 0, &numHexagons)); + H3Index *hexagonsTMH = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &transMeridianHoleGeoPolygon, 7, 0, hexagonsTMH)); + actualNumIndexes = countNonNullIndexes(hexagonsTMH, numHexagons); + + t_assert(actualNumIndexes == expectedSize - actualNumHoleIndexes, + "got expected polygonToCells size (transmeridian hole)"); + + free(hexagons); + free(hexagonsTM); + free(hexagonsTMFH); + free(hexagonsTMH); + } + + TEST(polygonToCellsTransmeridianComplex) { + // 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 + LatLng 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}}; + GeoLoop geoloop = {.numVerts = 6, .verts = verts}; + GeoPolygon polygon = {.geoloop = geoloop, .numHoles = 0}; + + int64_t numHexagons; + t_assertSuccess( + H3_EXPORT(maxPolygonToCellsSize)(&polygon, 4, 0, &numHexagons)); + + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + t_assertSuccess( + H3_EXPORT(polygonToCellsExperimental)(&polygon, 4, 0, hexagons)); + + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 1204, + "got expected polygonToCells size (complex transmeridian)"); + + free(hexagons); + } + + TEST(polygonToCellsPentagon) { + H3Index pentagon; + setH3Index(&pentagon, 9, 24, 0); + LatLng coord; + H3_EXPORT(cellToLatLng)(pentagon, &coord); + + // Length of half an edge of the polygon, in radians + double edgeLength2 = H3_EXPORT(degsToRads)(0.001); + + LatLng boundingTopRigt = coord; + boundingTopRigt.lat += edgeLength2; + boundingTopRigt.lng += edgeLength2; + + LatLng boundingTopLeft = coord; + boundingTopLeft.lat += edgeLength2; + boundingTopLeft.lng -= edgeLength2; + + LatLng boundingBottomRight = coord; + boundingBottomRight.lat -= edgeLength2; + boundingBottomRight.lng += edgeLength2; + + LatLng boundingBottomLeft = coord; + boundingBottomLeft.lat -= edgeLength2; + boundingBottomLeft.lng -= edgeLength2; + + LatLng verts[] = {boundingBottomLeft, boundingTopLeft, boundingTopRigt, + boundingBottomRight}; + + GeoLoop geoloop; + geoloop.verts = verts; + geoloop.numVerts = 4; + + GeoPolygon polygon; + polygon.geoloop = geoloop; + polygon.numHoles = 0; + + int64_t numHexagons; + t_assertSuccess( + H3_EXPORT(maxPolygonToCellsSize)(&polygon, 9, 0, &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess( + H3_EXPORT(polygonToCellsExperimental)(&polygon, 9, 0, hexagons)); + + int found = 0; + int numPentagons = 0; + for (int i = 0; i < numHexagons; i++) { + if (hexagons[i] != 0) { + found++; + } + if (H3_EXPORT(isPentagon)(hexagons[i])) { + numPentagons++; + } + } + t_assert(found == 1, "one index found"); + t_assert(numPentagons == 1, "one pentagon found"); + free(hexagons); + } + + TEST(invalidFlags) { + int64_t numHexagons; + for (uint32_t flags = 1; flags <= 32; flags++) { + t_assert( + H3_EXPORT(maxPolygonToCellsSize)( + &sfGeoPolygon, 9, flags, &numHexagons) == E_OPTION_INVALID, + "Flags other than 0 are invalid for maxPolygonToCellsSize"); + } + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + for (uint32_t flags = 1; flags <= 32; flags++) { + t_assert(H3_EXPORT(polygonToCellsExperimental)( + &sfGeoPolygon, 9, flags, hexagons) == E_OPTION_INVALID, + "Flags other than 0 are invalid for polygonToCells"); + } + free(hexagons); + } + + TEST(fillIndex) { + iterateAllIndexesAtRes(0, fillIndex_assertions); + iterateAllIndexesAtRes(1, fillIndex_assertions); + iterateAllIndexesAtRes(2, fillIndex_assertions); + } +} diff --git a/src/apps/testapps/testPolygonToCellsReportedExperimental.c b/src/apps/testapps/testPolygonToCellsReportedExperimental.c new file mode 100644 index 000000000..58e688504 --- /dev/null +++ b/src/apps/testapps/testPolygonToCellsReportedExperimental.c @@ -0,0 +1,207 @@ +/* + * Copyright 2017-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include + +#include "algos.h" +#include "constants.h" +#include "h3Index.h" +#include "latLng.h" +#include "polyfill.h" +#include "test.h" +#include "utility.h" + +// Tests for specific polygonToCells examples + +SUITE(polygonToCells_reported) { + // https://github.com/uber/h3-js/issues/76#issuecomment-561204505 + TEST(entireWorld) { + // TODO: Fails for a single worldwide polygon + LatLng worldVerts[] = { + {-M_PI_2, -M_PI}, {M_PI_2, -M_PI}, {M_PI_2, 0}, {-M_PI_2, 0}}; + GeoLoop worldGeoLoop = {.numVerts = 4, .verts = worldVerts}; + GeoPolygon worldGeoPolygon = {.geoloop = worldGeoLoop, .numHoles = 0}; + LatLng worldVerts2[] = { + {-M_PI_2, 0}, {M_PI_2, 0}, {M_PI_2, M_PI}, {-M_PI_2, M_PI}}; + GeoLoop worldGeoLoop2 = {.numVerts = 4, .verts = worldVerts2}; + GeoPolygon worldGeoPolygon2 = {.geoloop = worldGeoLoop2, .numHoles = 0}; + + for (int res = 0; res < 3; res++) { + int64_t polygonToCellsSize; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)( + &worldGeoPolygon, res, 0, &polygonToCellsSize)); + H3Index *polygonToCellsOut = + calloc(polygonToCellsSize, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &worldGeoPolygon, res, 0, polygonToCellsOut)); + int64_t actualNumIndexes = + countNonNullIndexes(polygonToCellsOut, polygonToCellsSize); + + int64_t polygonToCellsSize2; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)( + &worldGeoPolygon2, res, 0, &polygonToCellsSize2)); + H3Index *polygonToCellsOut2 = + calloc(polygonToCellsSize2, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &worldGeoPolygon2, res, 0, polygonToCellsOut2)); + int64_t actualNumIndexes2 = + countNonNullIndexes(polygonToCellsOut2, polygonToCellsSize2); + + int64_t expectedTotalWorld; + t_assertSuccess(H3_EXPORT(getNumCells)(res, &expectedTotalWorld)); + t_assert(actualNumIndexes + actualNumIndexes2 == expectedTotalWorld, + "got expected polygonToCells size (entire world)"); + + // Sets should be disjoint + for (int i = 0; i < polygonToCellsSize; i++) { + if (polygonToCellsOut[i] == 0) continue; + + bool found = false; + for (int j = 0; j < polygonToCellsSize2; j++) { + if (polygonToCellsOut[i] == polygonToCellsOut2[j]) { + found = true; + break; + } + } + t_assert( + !found, + "Index found more than once when polygonToCellsing the " + "entire world"); + } + + free(polygonToCellsOut); + free(polygonToCellsOut2); + } + } + + // https://github.com/uber/h3-js/issues/67 + TEST(h3js_67) { + double east = H3_EXPORT(degsToRads)(-56.25); + double north = H3_EXPORT(degsToRads)(-33.13755119234615); + double south = H3_EXPORT(degsToRads)(-34.30714385628804); + double west = H3_EXPORT(degsToRads)(-57.65625); + + LatLng testVerts[] = { + {north, east}, {south, east}, {south, west}, {north, west}}; + GeoLoop testGeoLoop = {.numVerts = 4, .verts = testVerts}; + GeoPolygon testPolygon; + testPolygon.geoloop = testGeoLoop; + testPolygon.numHoles = 0; + + int res = 7; + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res, + 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 4499, + "got expected polygonToCells size (h3-js#67)"); + free(hexagons); + } + + // 2nd test case from h3-js#67 + TEST(h3js_67_2nd) { + double east = H3_EXPORT(degsToRads)(-57.65625); + double north = H3_EXPORT(degsToRads)(-34.30714385628804); + double south = H3_EXPORT(degsToRads)(-35.4606699514953); + double west = H3_EXPORT(degsToRads)(-59.0625); + + LatLng testVerts[] = { + {north, east}, {south, east}, {south, west}, {north, west}}; + GeoLoop testGeoLoop = {.numVerts = 4, .verts = testVerts}; + GeoPolygon testPolygon; + testPolygon.geoloop = testGeoLoop; + testPolygon.numHoles = 0; + + int res = 7; + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res, + 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + t_assert(actualNumIndexes == 4609, + "got expected polygonToCells size (h3-js#67, 2nd case)"); + free(hexagons); + } + + // https://github.com/uber/h3/issues/136 + TEST(h3_136) { + LatLng testVerts[] = {{0.10068990369902957, 0.8920772174196191}, + {0.10032914690616246, 0.8915914753447348}, + {0.10033349237998787, 0.8915860128746426}, + {0.10069496685903621, 0.8920742194546231}}; + GeoLoop testGeoLoop = {.numVerts = 4, .verts = testVerts}; + GeoPolygon testPolygon; + testPolygon.geoloop = testGeoLoop; + testPolygon.numHoles = 0; + + int res = 13; + int64_t numHexagons; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0, + &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res, + 0, hexagons)); + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + + 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(polygonToCellsExperimental)(&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/bbox.h b/src/h3lib/include/bbox.h index e81912b19..3ae68cbf8 100644 --- a/src/h3lib/include/bbox.h +++ b/src/h3lib/include/bbox.h @@ -22,6 +22,7 @@ #include +#include "h3api.h" #include "latLng.h" /** @struct BBox @@ -34,12 +35,20 @@ typedef struct { double west; ///< west longitude } BBox; +double bboxWidthRads(const BBox *bbox); +double bboxHeightRads(const BBox *bbox); bool bboxIsTransmeridian(const BBox *bbox); void bboxCenter(const BBox *bbox, LatLng *center); bool bboxContains(const BBox *bbox, const LatLng *point); +bool bboxContainsBBox(const BBox *a, const BBox *b); +bool bboxOverlapsBBox(const BBox *a, const BBox *b); bool bboxEquals(const BBox *b1, const BBox *b2); H3Error bboxHexEstimate(const BBox *bbox, int res, int64_t *out); H3Error lineHexEstimate(const LatLng *origin, const LatLng *destination, int res, int64_t *out); +void scaleBBox(BBox *bbox, double scale); +void bboxNormalization(const BBox *a, const BBox *b, + LongitudeNormalization *aNormalization, + LongitudeNormalization *bNormalization); #endif diff --git a/src/h3lib/include/iterators.h b/src/h3lib/include/iterators.h index ad735ed6f..9e1de9bfd 100644 --- a/src/h3lib/include/iterators.h +++ b/src/h3lib/include/iterators.h @@ -54,6 +54,7 @@ typedef struct { } IterCellsChildren; DECLSPEC IterCellsChildren iterInitParent(H3Index h, int childRes); +void _iterInitParent(H3Index h, int childRes, IterCellsChildren *iter); DECLSPEC IterCellsChildren iterInitBaseCellNum(int baseCellNum, int childRes); DECLSPEC void iterStepChild(IterCellsChildren *iter); diff --git a/src/h3lib/include/latLng.h b/src/h3lib/include/latLng.h index 2dc674efa..8d644e578 100644 --- a/src/h3lib/include/latLng.h +++ b/src/h3lib/include/latLng.h @@ -32,9 +32,17 @@ /** epsilon of ~0.1mm in radians */ #define EPSILON_RAD (EPSILON_DEG * M_PI_180) +typedef enum { + NORMALIZE_NONE = 0, // Do not normalize + NORMALIZE_EAST = 1, // Normalize negative numbers to the east + NORMALIZE_WEST = 2 // Normalize positive numbers to the west +} LongitudeNormalization; + void setGeoDegs(LatLng *p, double latDegs, double lngDegs); double constrainLat(double lat); double constrainLng(double lng); +double normalizeLng(const double lng, + const LongitudeNormalization normalization); bool geoAlmostEqual(const LatLng *p1, const LatLng *p2); bool geoAlmostEqualThreshold(const LatLng *p1, const LatLng *p2, diff --git a/src/h3lib/include/polyfill.h b/src/h3lib/include/polyfill.h new file mode 100644 index 000000000..3bed9bc67 --- /dev/null +++ b/src/h3lib/include/polyfill.h @@ -0,0 +1,87 @@ +/* + * Copyright 2023 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file polyfill.h + * @brief Functions relating to the polygon-to-cells functionality + */ + +#ifndef POLYFILL_H +#define POLYFILL_H + +#include + +#include "bbox.h" +#include "h3api.h" +#include "iterators.h" + +/** + * IterCellsPolygonCompact: struct for iterating through all the cells + * within a given polygon, outputting a compact set. + * + * Constructors: + * + * Initialize with `iterInitPolygon`. This will save a reference to the + * input polygon and allocate memory for data structures used in the + * iteration. Iterators initialized in this way must be destroyed by + * `iterDestroyPolygon` to free allocated memory. + * + * Iteration: + * + * Step iterator with `iterStepPolygonCompact`. + * During the lifetime of the `iterCellsPolygonCompact`, the current iterate + * is accessed via the `iterCellsPolygonCompact.cell` member. When the iterator + * is exhausted or if there was an error in initialization or iteration, + * `iterCellsPolygonCompact.cell` will be `H3_NULL` after calling + * `iterStepChild`. It is the responsibiliy of the caller to check + * `iterCellsPolygonCompact.error` when `H3_NULL` is received. + * + * Cleanup: + * + * Destroy the iterator and free allocated memory with `iterDestroyPolygon`. + */ +typedef struct { + H3Index cell; // current value + H3Error error; // error, if any + int _res; // target resolution + uint32_t _flags; // Mode flags for the polygonToCells operation + const GeoPolygon *_polygon; // the polygon we're filling + BBox *_bboxes; // Bounding box(es) for the polygon and its holes + bool _started; // Whether iteration has started +} IterCellsPolygonCompact; + +DECLSPEC IterCellsPolygonCompact +iterInitPolygonCompact(const GeoPolygon *polygon, int res, uint32_t flags); +DECLSPEC void iterStepPolygonCompact(IterCellsPolygonCompact *iter); +DECLSPEC void iterDestroyPolygonCompact(IterCellsPolygonCompact *iter); + +typedef struct { + H3Index cell; // current value + H3Error error; // error, if any + IterCellsPolygonCompact _cellIter; // sub-iterator for compact cells + IterCellsChildren _childIter; // sub-iterator for cell children +} IterCellsPolygon; + +DECLSPEC IterCellsPolygon iterInitPolygon(const GeoPolygon *polygon, int res, + uint32_t flags); +DECLSPEC void iterStepPolygon(IterCellsPolygon *iter); +DECLSPEC void iterDestroyPolygon(IterCellsPolygon *iter); + +DECLSPEC H3Error H3_EXPORT(polygonToCellsExperimental)( + const GeoPolygon *polygon, int res, uint32_t flags, H3Index *out); + +H3Error cellToBBox(H3Index cell, BBox *out, bool coverChildren); +H3Index baseCellNumToCell(int baseCellNum); + +#endif diff --git a/src/h3lib/include/polygon.h b/src/h3lib/include/polygon.h index 3641f9fbe..387ed034a 100644 --- a/src/h3lib/include/polygon.h +++ b/src/h3lib/include/polygon.h @@ -44,6 +44,14 @@ void bboxesFromGeoPolygon(const GeoPolygon *polygon, BBox *bboxes); bool pointInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes, const LatLng *coord); +bool cellBoundaryInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes, + const CellBoundary *boundary, + const BBox *boundaryBBox); +bool cellBoundaryCrossesGeoLoop(const GeoLoop *geoloop, const BBox *loopBBox, + const CellBoundary *boundary, + const BBox *boundaryBBox); +bool lineIntersectsLine(const LatLng *a1, const LatLng *a2, const LatLng *b1, + const LatLng *b2); // The following functions are created via macro in polygonAlgos.h, // so their signatures are documented here: diff --git a/src/h3lib/lib/bbox.c b/src/h3lib/lib/bbox.c index 033b02edc..b86d88000 100644 --- a/src/h3lib/lib/bbox.c +++ b/src/h3lib/lib/bbox.c @@ -27,6 +27,23 @@ #include "h3Index.h" #include "latLng.h" +/** + * Width of the bounding box, in rads + * @param bbox Bounding box to inspect + * @return width, in rads + */ +double bboxWidthRads(const BBox *bbox) { + return bboxIsTransmeridian(bbox) ? bbox->east - bbox->west + M_2PI + : bbox->east - bbox->west; +} + +/** + * Height of the bounding box + * @param bbox Bounding box to inspect + * @return height, in rads + */ +double bboxHeightRads(const BBox *bbox) { return bbox->north - bbox->south; } + /** * Whether the given bounding box crosses the antimeridian * @param bbox Bounding box to inspect @@ -62,6 +79,56 @@ bool bboxContains(const BBox *bbox, const LatLng *point) { (point->lng >= bbox->west && point->lng <= bbox->east)); } +/** + * Whether two bounding boxes overlap + * @param a First bounding box + * @param b Second bounding box + * @return Whether the bounding boxes overlap + */ +bool bboxOverlapsBBox(const BBox *a, const BBox *b) { + // Check whether latitude coords overlap + if (a->north < b->south || a->south > b->north) { + return false; + } + + // Check whether longitude coords overlap, accounting for transmeridian + // bboxes + LongitudeNormalization aNormalization; + LongitudeNormalization bNormalization; + bboxNormalization(a, b, &aNormalization, &bNormalization); + + if (normalizeLng(a->east, aNormalization) < + normalizeLng(b->west, bNormalization) || + normalizeLng(a->west, aNormalization) > + normalizeLng(b->east, bNormalization)) { + return false; + } + + return true; +} + +/** + * Whether one bounding box contains another + * @param a First bounding box + * @param b Second bounding box + * @return Whether a contains b + */ +bool bboxContainsBBox(const BBox *a, const BBox *b) { + // Check whether latitude coords are contained + if (a->north < b->north || a->south > b->south) { + return false; + } + // Check whether longitude coords are contained + // Account for transmeridian bboxes + LongitudeNormalization aNormalization; + LongitudeNormalization bNormalization; + bboxNormalization(a, b, &aNormalization, &bNormalization); + return normalizeLng(a->west, aNormalization) <= + normalizeLng(b->west, bNormalization) && + normalizeLng(a->east, aNormalization) >= + normalizeLng(b->east, bNormalization); +} + /** * Whether two bounding boxes are strictly equal * @param b1 Bounding box 1 @@ -175,3 +242,60 @@ H3Error lineHexEstimate(const LatLng *origin, const LatLng *destination, *out = estimate; return E_SUCCESS; } + +/** + * Scale a given bounding box by some factor. Scales both width and height + * by the factor, rather than scaling area, which will scale at scale^2. + * Note that this function is meant to handle bounding boxes and scales, + * within a reasonable domain, and does not guarantee reasonable results for + * extreme values. + * @param bbox Bounding box to scale, in-place + * @param scale Scale factor + */ +void scaleBBox(BBox *bbox, double scale) { + double width = bboxWidthRads(bbox); + double height = bboxHeightRads(bbox); + double widthBuffer = (width * scale - width) / 2; + double heightBuffer = (height * scale - height) / 2; + // Scale north and south, clamping to latitude domain + bbox->north += heightBuffer; + if (bbox->north > M_PI_2) bbox->north = M_PI_2; + bbox->south -= heightBuffer; + if (bbox->south < -M_PI_2) bbox->south = -M_PI_2; + // Scale east and west, clamping to longitude domain + bbox->east += widthBuffer; + if (bbox->east > M_PI) bbox->east -= M_2PI; + if (bbox->east < -M_PI) bbox->east += M_2PI; + bbox->west -= widthBuffer; + if (bbox->west > M_PI) bbox->west -= M_2PI; + if (bbox->west < -M_PI) bbox->west += M_2PI; +} + +/** + * Determine the longitude normalization scheme for two bounding boxes, either + * or both of which might cross the antimeridian. The goal is to transform + * latitudes in one or both boxes so that they are in the same frame of + * reference and can be operated on with standard Cartesian functions. + * @param a First bounding box + * @param b Second bounding box + * @param aNormalization Output: Normalization for longitudes in the first box + * @param bNormalization Output: Normalization for longitudes in the second box + */ +void bboxNormalization(const BBox *a, const BBox *b, + LongitudeNormalization *aNormalization, + LongitudeNormalization *bNormalization) { + bool aIsTransmeridian = bboxIsTransmeridian(a); + bool bIsTransmeridian = bboxIsTransmeridian(b); + bool aToBTrendsEast = a->west - b->east < b->west - a->east; + // If neither is transmeridian, no normalization. + // If both are transmeridian, normalize east by convention. + // If one is transmeridian and one is not, normalize toward the other. + *aNormalization = !aIsTransmeridian ? NORMALIZE_NONE + : bIsTransmeridian ? NORMALIZE_EAST + : aToBTrendsEast ? NORMALIZE_EAST + : NORMALIZE_WEST; + *bNormalization = !bIsTransmeridian ? NORMALIZE_NONE + : aIsTransmeridian ? NORMALIZE_EAST + : aToBTrendsEast ? NORMALIZE_WEST + : NORMALIZE_EAST; +} \ No newline at end of file diff --git a/src/h3lib/lib/iterators.c b/src/h3lib/lib/iterators.c index ac68fd1ef..ef4d1a785 100644 --- a/src/h3lib/lib/iterators.c +++ b/src/h3lib/lib/iterators.c @@ -213,27 +213,33 @@ and so on. */ IterCellsChildren iterInitParent(H3Index h, int childRes) { IterCellsChildren it; + _iterInitParent(h, childRes, &it); + return it; +} - it._parentRes = H3_GET_RESOLUTION(h); +/** + * Internal function - initialize a parent iterator in-place + */ +void _iterInitParent(H3Index h, int childRes, IterCellsChildren *iter) { + iter->_parentRes = H3_GET_RESOLUTION(h); - if (childRes < it._parentRes || childRes > MAX_H3_RES || h == H3_NULL) { - return _null_iter(); + if (childRes < iter->_parentRes || childRes > MAX_H3_RES || h == H3_NULL) { + *iter = _null_iter(); + return; } - it.h = _zeroIndexDigits(h, it._parentRes + 1, childRes); - H3_SET_RESOLUTION(it.h, childRes); + iter->h = _zeroIndexDigits(h, iter->_parentRes + 1, childRes); + H3_SET_RESOLUTION(iter->h, childRes); - if (H3_EXPORT(isPentagon)(it.h)) { + if (H3_EXPORT(isPentagon)(iter->h)) { // The skip digit skips `1` for pentagons. // The "_skipDigit" moves to the left as we count up from the // child resolution to the parent resolution. - it._skipDigit = childRes; + iter->_skipDigit = childRes; } else { // if not a pentagon, we can ignore "skip digit" logic - it._skipDigit = -1; + iter->_skipDigit = -1; } - - return it; } /** diff --git a/src/h3lib/lib/latLng.c b/src/h3lib/lib/latLng.c index 5c2c098e0..8d8e8cad2 100644 --- a/src/h3lib/lib/latLng.c +++ b/src/h3lib/lib/latLng.c @@ -137,6 +137,24 @@ double constrainLng(double lng) { return lng; } +/** + * Normalize an input longitude according to the specified normalization + * @param lng Input longitude + * @param normalization Longitude normalization strategy + * @return Normalized longitude + */ +double normalizeLng(const double lng, + const LongitudeNormalization normalization) { + switch (normalization) { + case NORMALIZE_EAST: + return lng < 0 ? lng + (double)M_2PI : lng; + case NORMALIZE_WEST: + return lng > 0 ? lng - (double)M_2PI : lng; + default: + return lng; + } +} + /** * The great circle distance in radians between two spherical coordinates. * diff --git a/src/h3lib/lib/polyfill.c b/src/h3lib/lib/polyfill.c new file mode 100644 index 000000000..7ba5a28b5 --- /dev/null +++ b/src/h3lib/lib/polyfill.c @@ -0,0 +1,593 @@ +/* + * Copyright 2023 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file polyfill.c + * @brief Functions relating to the cell-to-polygon algorithm + */ + +#include "polyfill.h" + +#include +#include + +#include "alloc.h" +#include "baseCells.h" +#include "coordijk.h" +#include "h3Assert.h" +#include "h3Index.h" +#include "polygon.h" + +// Factor by which to scale the cell bounding box to include all cells. +// This was determined empirically by finding the smallest factor that +// passed exhaustive tests. +#define CELL_SCALE_FACTOR 1.1 + +// Factor by which to scale the cell bounding box to include all children. +// This was determined empirically by finding the smallest factor that +// passed exhaustive tests. +#define CHILD_SCALE_FACTOR 1.4 + +/** + * Max cell edge length, in radians, for each resolution. This was computed + * by taking the max exact edge length for cells at the center of each base + * cell at that resolution. + */ +static double MAX_EDGE_LENGTH_RADS[MAX_H3_RES + 1] = { + 0.21577206265130, 0.08308767068495, 0.03148970436439, 0.01190662871439, + 0.00450053330908, 0.00170105523619, 0.00064293917678, 0.00024300820659, + 0.00009184847087, 0.00003471545901, 0.00001312121017, 0.00000495935129, + 0.00000187445860, 0.00000070847876, 0.00000026777980, 0.00000010121125}; + +/** All cells that contain the north pole, by res */ +static H3Index NORTH_POLE_CELLS[MAX_H3_RES + 1] = { + 0x8001fffffffffff, 0x81033ffffffffff, 0x820327fffffffff, 0x830326fffffffff, + 0x8403263ffffffff, 0x85032623fffffff, 0x860326237ffffff, 0x870326233ffffff, + 0x880326233bfffff, 0x890326233abffff, 0x8a0326233ab7fff, 0x8b0326233ab0fff, + 0x8c0326233ab03ff, 0x8d0326233ab03bf, 0x8e0326233ab039f, 0x8f0326233ab0399}; + +/** All cells that contain the south pole, by res */ +static H3Index SOUTH_POLE_CELLS[MAX_H3_RES + 1] = { + 0x80f3fffffffffff, 0x81f2bffffffffff, 0x82f297fffffffff, 0x83f293fffffffff, + 0x84f2939ffffffff, 0x85f29383fffffff, 0x86f29380fffffff, 0x87f29380effffff, + 0x88f29380e1fffff, 0x89f29380e0fffff, 0x8af29380e0d7fff, 0x8bf29380e0d0fff, + 0x8cf29380e0d0dff, 0x8df29380e0d0cff, 0x8ef29380e0d0cc7, 0x8ff29380e0d0cc4}; + +/** Pre-calculated bounding boxes for all res 0 cells */ +static BBox RES0_BBOXES[NUM_BASE_CELLS] = { + {1.52480158339146, 1.20305471830087, -0.60664883654036, 0.00568297271999}, + {1.52480158339146, 1.17872424267511, -0.60664883654036, 2.54046980298264}, + {1.52480158339146, 1.09069387298096, -2.85286053297673, 1.64310689027893}, + {1.41845302535151, 1.01285145697208, 0.00568297271999, -1.16770379632602}, + {1.27950477868453, 0.97226652536306, 0.55556064983494, -0.18229924845326}, + {1.32929586572429, 0.91898920750071, 2.05622344943192, 1.08813154278274}, + {1.32899086063916, 0.94271815376360, -2.29875289606378, 3.01700008041993}, + {1.26020983864103, 0.84291228415618, -0.89971867664861, -1.75967359310997}, + {1.21114673854945, 0.86170600921069, 1.19129757609455, 0.43777608996454}, + {1.21075831414294, 0.83795331049498, -1.72022875779891, -2.43793861727138}, + {1.15546530929588, 0.78982455384253, 2.53659412229266, 1.85709133451243}, + {1.15528445067052, 0.76641428724335, -3.06738507202411, 2.53646110244042}, + {1.10121643537669, 0.71330093663066, 0.09640581900154, -0.52154514518248}, + {1.07042472765165, 0.67603948819406, -0.47984202840088, -1.10306159603090}, + {1.03270228748960, 0.72356358827215, -2.24990138725146, -2.74510220919157}, + {1.01929924623886, 0.65491232835426, 0.63035574240731, 0.03537030096470}, + {1.01786037568858, 0.58827636737638, 1.53192721817065, 0.93672682511233}, + {0.98081434136020, 0.61076063532947, -2.67100636598529, 3.06516463008733}, + {0.98106023192774, 0.58679836571570, 2.02829766214461, 1.51334374970280}, + {0.96374551790056, 0.55186491737474, -1.42976721313659, -1.96852202530104}, + {0.87536136210723, 0.50008952762292, -1.92435613571430, -2.41641343219793}, + {0.88611243445554, 0.52742963716774, -0.95781946324194, -1.47628966305930}, + {0.86881343251986, 0.50770567021439, 1.03236795495839, 0.50347284027426}, + {0.89235638181782, 0.48781264892508, 2.76430302119150, 2.29989716697031}, + {0.82570569254601, 0.52173101741059, 2.30921681461428, 1.93198541828980}, + {0.80599330438546, 0.40150819579319, -3.06417559403240, 2.70079300784409}, + {0.81612079704781, 0.38396800633226, -0.21614378891839, -0.70420149722178}, + {0.75822779851431, 0.39943555383751, -2.34059978084699, -2.82127373822444}, + {0.78861390967531, 0.38742018303868, 0.23115687731652, -0.22599491086066}, + {0.71515840341957, 0.33012478438475, -0.64847976163163, -1.08249728121219}, + {0.70359051048414, 0.29148673180722, 1.71441081857246, 1.28443348381696}, + {0.69190629544818, 0.28808313184381, 0.64863909244647, 0.16372369282557}, + {0.64863235654749, 0.26290420067147, 2.10318098268379, 1.69556122548344}, + {0.65722892279906, 0.28222653310929, 1.30918693285466, 0.87594416271685}, + {0.64750997738584, 0.24149865709850, -1.30272192474556, -1.68708570163242}, + {0.62380174028378, 0.25522080363509, -2.72428423026826, 3.10401473237630}, + {0.64228460410023, 0.21206753429148, -1.67639240992071, -2.11772366767341}, + {0.59919175361146, 0.21620460836570, 2.48592868387690, 2.07350353893591}, + {0.55637406851384, 0.25276557437230, -0.99885388505694, -1.32642489358939}, + {0.55648013300665, 0.15187401321019, 2.87032088421324, 2.44642320475367}, + {0.54603687970450, 0.15589091511369, -2.06789866067060, -2.49091419631961}, + {0.51206347752697, 0.15522020377124, 0.95446767315996, 0.54443262110414}, + {0.49767951537101, 0.10944898890579, -0.04335162263358, -0.42900268178569}, + {0.46538045483671, 0.06029968637720, -0.41240613713421, -0.80603623808166}, + {0.44686891066946, 0.06926857458503, 0.32053284794952, -0.07005748900849}, + {0.43208958202064, 0.07796440938140, -3.06232453079660, 2.80602499990282}, + {0.43103892586713, 0.02927431919853, -2.41589238618422, -2.85735809951951}, + {0.38073727558986, -0.00297016159959, -0.77039553861218, -1.14788248745028}, + {0.39113816687141, -0.01518764903038, 1.49130246958290, 1.14714731736311}, + {0.33421063142418, 0.02526613430348, 1.15141032578749, 0.85000706261644}, + {0.38915669778582, -0.04371359825454, 1.88046353933242, 1.48230231380717}, + {0.33787520825987, -0.04835090128296, -1.12274014380603, -1.49454408844749}, + {0.33601418932337, -0.06675068178541, 2.23792354204464, 1.85723423013211}, + {0.31838318078049, -0.05821955623722, 0.66058854060373, 0.25452572938783}, + {0.33630761471457, -0.07589541031521, -1.47957331741818, -1.85981735718264}, + {0.28924817322870, -0.09150638064667, -1.83561930288569, -2.21855897384292}, + {0.26678632252475, -0.10058088990867, -2.76808651991421, 3.12792953247061}, + {0.29285254112587, -0.13483165093783, 2.61406468380434, 2.20466422911705}, + {0.20150342788824, -0.10279852729762, 0.06881896344365, -0.23925229432978}, + {0.21283813275258, -0.18626835417891, 2.93800440256577, 2.57470747655623}, + {0.19587614179884, -0.17237030304155, -2.16941795427335, -2.55405165906601}, + {0.17237030304155, -0.19587614179884, 0.97217469931645, 0.58754099452378}, + {0.18626835417891, -0.21283813275258, -0.20358825102402, -0.56688517703356}, + {0.10279852729762, -0.20150342788824, -3.07277369014614, 2.90234035926002}, + {0.13483165093783, -0.29285254112587, -0.52752796978545, -0.93692842447275}, + {0.10058088990867, -0.26678632252475, 0.37350613367558, -0.01366312111919}, + {0.09150638064667, -0.28924817322870, 1.30597335070410, 0.92303367974687}, + {0.07589541031521, -0.33630761471457, 1.66201933617161, 1.28177529640715}, + {0.05821955623722, -0.31838318078049, -2.48100411298606, -2.88706692420196}, + {0.06675068178541, -0.33601418932337, -0.90366911154516, -1.28435842345769}, + {0.04835090128296, -0.33787520825987, 2.01885250978376, 1.64704856514230}, + {0.04371359825454, -0.38915669778582, -1.26112911425737, -1.65929033978262}, + {-0.02526613430348, -0.33421063142418, -1.99018232780231, + -2.29158559097336}, + {0.01518764903038, -0.39113816687140, -1.65029018400690, -1.99444533622668}, + {0.00297016159959, -0.38073727558986, 2.37119711497761, 1.99371016613951}, + {-0.02927431919853, -0.43103892586713, 0.72570026740558, 0.28423455407029}, + {-0.07796440938140, -0.43208958202064, 0.07926812279319, -0.33556765368697}, + {-0.06926857458503, -0.44686891066946, -2.82105980564027, 3.07153516458131}, + {-0.06029968637720, -0.46538045483671, 2.72918651645558, 2.33555641550814}, + {-0.10944898890579, -0.49767951537101, 3.09824103095621, 2.71258997180410}, + {-0.15522020377124, -0.51206347752697, -2.18712498042983, + -2.59716003248565}, + {-0.15589091511369, -0.54603687970450, 1.07369399291919, 0.65067845727018}, + {-0.15187401321019, -0.55648013300665, -0.27127176937655, + -0.69516944883612}, + {-0.25276557437230, -0.55637406851385, 2.14273876853285, 1.81516776000041}, + {-0.21620460836570, -0.59919175361146, -0.65566396971290, + -1.06808911465388}, + {-0.21206753429148, -0.64228460410023, 1.46520024366909, 1.02386898591638}, + {-0.25522080363509, -0.62380174028378, 0.41730842332153, -0.03757792121350}, + {-0.24149865709850, -0.64750997738584, 1.83887072884423, 1.45450695195737}, + {-0.28222653310929, -0.65722892279906, -1.83240572073513, + -2.26564849087294}, + {-0.26290420067147, -0.64863235654749, -1.03841167090601, + -1.44603142810635}, + {-0.28808313184381, -0.69190629544818, -2.49295356114332, + -2.97786896076422}, + {-0.29148673180722, -0.70359051048414, -1.42718183501734, + -1.85715916977284}, + {-0.33012478438475, -0.71515840341957, 2.49311289195816, 2.05909537237761}, + {-0.38742018303868, -0.78861390967531, -2.91043577627328, 2.91559774272914}, + {-0.39943555383751, -0.75822779851431, 0.80099287274280, 0.32031891536535}, + {-0.38396800633226, -0.81612079704781, 2.92544886467140, 2.43739115636801}, + {-0.40150819579319, -0.80599330438546, 0.07741705955739, -0.44079964574570}, + {-0.52173101741059, -0.82570569254601, -0.83237583897551, + -1.20960723529999}, + {-0.48781264892508, -0.89235638181782, -0.37728963239830, + -0.84169548661948}, + {-0.50770567021439, -0.86881343251986, -2.10922469863141, + -2.63811981331554}, + {-0.52742963716774, -0.88611243445554, 2.18377319034785, 1.66530299053050}, + {-0.50008952762292, -0.87536136210723, 1.21723651787549, 0.72517922139186}, + {-0.55186491737474, -0.96374551790056, 1.71182544045320, 1.17307062828876}, + {-0.58679836571570, -0.98106023192774, -1.11329499144518, + -1.62824890388699}, + {-0.61076063532947, -0.98081434136020, 0.47058628760450, -0.07642802350246}, + {-0.58827636737638, -1.01786037568858, -1.60966543541914, + -2.20486582847747}, + {-0.65491232835426, -1.01929924623886, -2.51123691118248, + -3.10622235262510}, + {-0.72356358827215, -1.03270228748960, 0.89169126633833, 0.39649044439822}, + {-0.67603948819406, -1.07042472765165, 2.66175062518892, 2.03853105755889}, + {-0.71330093663066, -1.10121643537669, -3.04518683458825, 2.62004750840731}, + {-0.76641428724335, -1.15528445067052, 0.07420758156568, -0.60513155114938}, + {-0.78982455384253, -1.15546530929588, -0.60499853129713, + -1.28450131907736}, + {-0.83795331049498, -1.21075831414294, 1.42136389579088, 0.70365403631841}, + {-0.86170600921069, -1.21114673854945, -1.95029507749525, + -2.70381656362525}, + {-0.84291228415618, -1.26020983864103, 2.24187397694118, 1.38191906047983}, + {-0.94271815376360, -1.32899086063916, 0.84283975752601, -0.12459257316986}, + {-0.91898920750071, -1.32929586572429, -1.08536920415787, + -2.05346111080706}, + {-0.97226652536306, -1.27950477868453, -2.58603200375485, 2.95929340513654}, + {-1.01285145697208, -1.41845302535151, -3.13590968086981, 1.97388885726377}, + {-1.09069387298096, -1.52480158339146, 0.28873212061306, -1.49848576331087}, + {-1.17872424267511, -1.52480158339146, 2.53494381704943, -0.60112285060716}, + {-1.20305471830087, -1.52480158339146, -0.60112285060716, + 2.53494381704943}}; + +/** + * For a given cell, return its bounding box. If coverChildren is true, the bbox + * will be guaranteed to contain its children at any finer resolution. Note that + * no guarantee is provided as to the level of accuracy, and the bounding box + * may have a significant margin of error. + * @param cell Cell to calculate bbox for + * @param out BBox to hold output + * @param coverChildren Whether the bounding box should cover all children + */ +H3Error cellToBBox(H3Index cell, BBox *out, bool coverChildren) { + // Adjust the BBox to handle poles, if needed + int res = H3_GET_RESOLUTION(cell); + + if (res == 0) { + int baseCell = H3_GET_BASE_CELL(cell); + if (NEVER(baseCell < 0) || baseCell >= NUM_BASE_CELLS) { + return E_CELL_INVALID; + } + *out = RES0_BBOXES[baseCell]; + } else { + LatLng center; + H3Error centerErr = H3_EXPORT(cellToLatLng)(cell, ¢er); + if (centerErr != E_SUCCESS) { + return centerErr; + } + double lngRatio = 1 / cos(center.lat); + out->north = center.lat + MAX_EDGE_LENGTH_RADS[res]; + out->south = center.lat - MAX_EDGE_LENGTH_RADS[res]; + out->east = center.lng + MAX_EDGE_LENGTH_RADS[res] * lngRatio; + out->west = center.lng - MAX_EDGE_LENGTH_RADS[res] * lngRatio; + } + + // Buffer the bounding box to cover children. Call this even if no buffering + // is required in order to normalize the bbox to lat/lng bounds + scaleBBox(out, coverChildren ? CHILD_SCALE_FACTOR : CELL_SCALE_FACTOR); + + // Cell that contains the north pole + if (cell == NORTH_POLE_CELLS[res]) { + out->north = M_PI_2; + } + + // Cell that contains the north pole + if (cell == SOUTH_POLE_CELLS[res]) { + out->south = -M_PI_2; + } + + // If we contain a pole, expand the longitude to include the full domain, + // effectively making the bbox a circle around the pole. + if (out->north == M_PI_2 || out->south == -M_PI_2) { + out->east = M_PI; + out->west = -M_PI; + } + + return E_SUCCESS; +} + +/** + * Get a base cell by number, or H3_NULL if out of bounds + */ +H3Index baseCellNumToCell(int baseCellNum) { + if (baseCellNum < 0 || baseCellNum >= NUM_BASE_CELLS) { + return H3_NULL; + } + H3Index baseCell; + setH3Index(&baseCell, 0, baseCellNum, 0); + return baseCell; +} + +static void iterErrorPolygonCompact(IterCellsPolygonCompact *iter, + H3Error error) { + iterDestroyPolygonCompact(iter); + iter->error = error; +} + +/** + * Given a cell, find the next cell in the sequence of all cells + * to check in the iteration. + */ +static H3Index nextCell(H3Index cell) { + int res = H3_GET_RESOLUTION(cell); + while (true) { + // If this is a base cell, set to next base cell (or H3_NULL if done) + if (res == 0) { + return baseCellNumToCell(H3_GET_BASE_CELL(cell) + 1); + } + + // Faster cellToParent when we know the resolution is valid + // and we're only moving up one level + H3Index parent = cell; + H3_SET_RESOLUTION(parent, res - 1); + H3_SET_INDEX_DIGIT(parent, res, H3_DIGIT_MASK); + + // If not the last sibling of parent, return next sibling + Direction digit = H3_GET_INDEX_DIGIT(cell, res); + if (digit < INVALID_DIGIT - 1) { + H3_SET_INDEX_DIGIT(cell, res, + digit + ((H3_EXPORT(isPentagon)(parent) && + digit == CENTER_DIGIT) + ? 2 // Skip missing pentagon child + : 1)); + return cell; + } + // Move up to the parent for the next loop iteration + res--; + cell = parent; + } +} + +/** + * Initialize a IterCellsPolygonCompact struct representing the sequence of + * compact cells within the target polygon. The test for including edge cells is + * defined by the polyfill mode passed in the `flags` argument. + * + * Initialization of this object may fail, in which case the `error` property + * will be set and all iteration will return H3_NULL. It is the responsibility + * of the caller to check the error property after initialization. + * + * At any point in the iteration, starting once the struct is initialized, the + * output value can be accessed through the `cell` property. + * + * Note that initializing the iterator allocates memory. If an iterator is + * exhausted or returns an error that memory is released; otherwise it must be + * released manually with iterDestroyPolygonCompact. + * + * @param polygon Polygon to fill with compact cells + * @param res Finest resolution for output cells + * @param flags Bit mask of option flags + * @return Initialized iterator, with the first value available + */ +IterCellsPolygonCompact iterInitPolygonCompact(const GeoPolygon *polygon, + int res, uint32_t flags) { + IterCellsPolygonCompact iter = {// Initialize output properties. The first + // valid cell will be set in iterStep + .cell = baseCellNumToCell(0), + .error = E_SUCCESS, + // Save input arguments + ._polygon = polygon, + ._res = res, + ._flags = flags, + ._bboxes = NULL, + ._started = false}; + + if (res < 0 || res > MAX_H3_RES) { + iterErrorPolygonCompact(&iter, E_RES_DOMAIN); + return iter; + } + + if (flags != 0) { + iterErrorPolygonCompact(&iter, E_OPTION_INVALID); + return iter; + } + + // Initialize bounding boxes for polygon and any holes. Memory allocated + // here must be released through iterDestroyPolygonCompact + iter._bboxes = H3_MEMORY(calloc)((polygon->numHoles + 1), sizeof(BBox)); + if (!iter._bboxes) { + iterErrorPolygonCompact(&iter, E_MEMORY_ALLOC); + return iter; + } + bboxesFromGeoPolygon(polygon, iter._bboxes); + + // Start the iterator by taking the first step. + // This is necessary to have a valid value after initialization. + iterStepPolygonCompact(&iter); + + return iter; +} + +/** + * Increment the polyfill iterator, running the polygon to cells algorithm. + * + * Briefly, the algorithm checks every cell in the global grid hierarchically, + * starting with the base cells. Cells coarser than the target resolution are + * checked for complete child inclusion using a bounding box guaranteed to + * contain all children. + * - If the bounding box is contained by the polygon, output is set to the cell + * - If the bounding box intersects, recurse into the first child + * - Otherwise, continue with the next cell in sequence + * + * For cells at the target resolution, a finer-grained check is used according + * to the inclusion criteria set in flags. + * + * @param iter Iterator to increment + */ +void iterStepPolygonCompact(IterCellsPolygonCompact *iter) { + H3Index cell = iter->cell; + + // once the cell is H3_NULL, the iterator returns an infinite sequence of + // H3_NULL + if (cell == H3_NULL) return; + + // For the first step, we need to evaluate the current cell; after that, we + // should start with the next cell. + if (iter->_started) { + cell = nextCell(cell); + } else { + iter->_started = true; + } + + while (cell) { + int cellRes = H3_GET_RESOLUTION(cell); + + // Target res: Do a fine-grained check + if (cellRes == iter->_res) { + // Check if the cell is in the polygon + // TODO: Handle other polyfill modes here + LatLng center; + H3Error centerErr = H3_EXPORT(cellToLatLng)(cell, ¢er); + if (NEVER(centerErr != E_SUCCESS)) { + iterErrorPolygonCompact(iter, centerErr); + return; + } + if (pointInsidePolygon(iter->_polygon, iter->_bboxes, ¢er)) { + // Set to next output + iter->cell = cell; + return; + } + } + + // Coarser cell: Check the bounding box + if (cellRes < iter->_res) { + // Get a bounding box for all of the cell's children + BBox bbox; + H3Error bboxErr = cellToBBox(cell, &bbox, true); + if (bboxErr) { + iterErrorPolygonCompact(iter, bboxErr); + return; + } + if (bboxOverlapsBBox(&iter->_bboxes[0], &bbox)) { + // Quick check for possible containment + if (bboxContainsBBox(&iter->_bboxes[0], &bbox)) { + // Convert bbox to cell boundary, CCW vertex order + CellBoundary bboxBoundary = { + .numVerts = 4, + .verts = {{bbox.north, bbox.east}, + {bbox.north, bbox.west}, + {bbox.south, bbox.west}, + {bbox.south, bbox.east}}}; + // Do a fine-grained, more expensive check on the polygon + if (cellBoundaryInsidePolygon(iter->_polygon, iter->_bboxes, + &bboxBoundary, &bbox)) { + // Bounding box is fully contained, so all children are + // included. Set to next output. + iter->cell = cell; + return; + } + } + // Otherwise, the intersecting bbox means we need to test all + // children, starting with the first child + H3Index child; + H3Error childErr = + H3_EXPORT(cellToCenterChild)(cell, cellRes + 1, &child); + if (childErr) { + iterErrorPolygonCompact(iter, childErr); + return; + } + // Restart the loop with the child cell + cell = child; + continue; + } + } + + // Find the next cell in the sequence of all cells and continue + cell = nextCell(cell); + } + // If we make it out of the loop, we're done + iterDestroyPolygonCompact(iter); +} + +/** + * Destroy an iterator, releasing any allocated memory. Iterators destroyed in + * this manner are safe to use but will always return H3_NULL. + * @param iter Iterator to destroy + */ +void iterDestroyPolygonCompact(IterCellsPolygonCompact *iter) { + if (iter->_bboxes) { + H3_MEMORY(free)(iter->_bboxes); + } + iter->cell = H3_NULL; + iter->error = E_SUCCESS; + iter->_polygon = NULL; + iter->_res = -1; + iter->_flags = 0; + iter->_bboxes = NULL; +} + +/** + * Initialize a IterCellsPolygon struct representing the sequence of + * cells within the target polygon. The test for including edge cells is defined + * by the polyfill mode passed in the `flags` argument. + * + * Initialization of this object may fail, in which case the `error` property + * will be set and all iteration will return H3_NULL. It is the responsibility + * of the caller to check the error property after initialization. + * + * At any point in the iteration, starting once the struct is initialized, the + * output value can be accessed through the `cell` property. + * + * Note that initializing the iterator allocates memory. If an iterator is + * exhausted or returns an error that memory is released; otherwise it must be + * released manually with iterDestroyPolygon. + * + * @param polygon Polygon to fill with cells + * @param res Resolution for output cells + * @param flags Bit mask of option flags + * @return Initialized iterator, with the first value available + */ +IterCellsPolygon iterInitPolygon(const GeoPolygon *polygon, int res, + uint32_t flags) { + // Create the sub-iterator for compact cells + IterCellsPolygonCompact cellIter = + iterInitPolygonCompact(polygon, res, flags); + // Create the sub-iterator for children + IterCellsChildren childIter = iterInitParent(cellIter.cell, res); + + IterCellsPolygon iter = {.cell = childIter.h, + .error = cellIter.error, + ._cellIter = cellIter, + ._childIter = childIter}; + return iter; +} + +/** + * Increment the polyfill iterator, outputting the latest cell at the + * desired resolution. + * + * @param iter Iterator to increment + */ +void iterStepPolygon(IterCellsPolygon *iter) { + if (iter->cell == H3_NULL) return; + + // See if there are more children to output + iterStepChild(&(iter->_childIter)); + if (iter->_childIter.h) { + iter->cell = iter->_childIter.h; + return; + } + + // Otherwise, increment the polyfill iterator + iterStepPolygonCompact(&(iter->_cellIter)); + if (iter->_cellIter.cell) { + _iterInitParent(iter->_cellIter.cell, iter->_cellIter._res, + &(iter->_childIter)); + iter->cell = iter->_childIter.h; + return; + } + + // All done, set to null and report errors if any + iter->cell = H3_NULL; + iter->error = iter->_cellIter.error; +} + +/** + * Destroy an iterator, releasing any allocated memory. Iterators destroyed in + * this manner are safe to use but will always return H3_NULL. + * @param iter Iterator to destroy + */ +void iterDestroyPolygon(IterCellsPolygon *iter) { + iterDestroyPolygonCompact(&(iter->_cellIter)); + // null out the child iterator by passing H3_NULL + _iterInitParent(H3_NULL, 0, &(iter->_childIter)); + iter->cell = H3_NULL; + iter->error = E_SUCCESS; +} + +/** + * polygonToCells takes a given GeoJSON-like data structure and preallocated, + * zeroed memory, and fills it with the hexagons that are contained by + * the GeoJSON-like data structure. Polygons are considered in Cartesian space. + * + * @param geoPolygon The geoloop and holes defining the relevant area + * @param res The Hexagon resolution (0-15) + * @param out The slab of zeroed memory to write to. Assumed to be big enough. + */ +H3Error H3_EXPORT(polygonToCellsExperimental)(const GeoPolygon *polygon, + int res, uint32_t flags, + H3Index *out) { + IterCellsPolygon iter = iterInitPolygon(polygon, res, flags); + int64_t i = 0; + for (; iter.cell; iterStepPolygon(&iter)) { + out[i++] = iter.cell; + } + return iter.error; +} \ No newline at end of file diff --git a/src/h3lib/lib/polygon.c b/src/h3lib/lib/polygon.c index ac880e7a3..a92fecfba 100644 --- a/src/h3lib/lib/polygon.c +++ b/src/h3lib/lib/polygon.c @@ -83,3 +83,124 @@ bool pointInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes, return contains; } + +/** + * Whether a cell boundary is completely contained by a polygon. + * @param geoPolygon The polygon to test + * @param bboxes The bboxes for the main geoloop and each of its holes + * @param boundary The cell boundary to test + * @return Whether the cell boundary is contained + */ +bool cellBoundaryInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes, + const CellBoundary *boundary, + const BBox *boundaryBBox) { + // First test a single point. Note that this fails fast if point is outside + // bounding box. + if (!pointInsidePolygon(geoPolygon, &bboxes[0], &boundary->verts[0])) { + return false; + } + + // If a point is contained, check for any line intersections + if (cellBoundaryCrossesGeoLoop(&(geoPolygon->geoloop), &bboxes[0], boundary, + boundaryBBox)) { + return false; + } + // Check for line intersections with any hole + for (int i = 0; i < geoPolygon->numHoles; i++) { + if (cellBoundaryCrossesGeoLoop(&(geoPolygon->holes[i]), &bboxes[i + 1], + boundary, boundaryBBox)) { + return false; + } + } + return true; +} + +/** + * Whether a cell boundary crosses a geo loop. Crossing in this case means + * whether any line segments intersect; it does not include containment. + * @param geoloop Geo loop to test + * @param boundary Cell boundary to test + * @return Whether any line segments in the boundary intersect any line + * segments in the geo loop + */ +bool cellBoundaryCrossesGeoLoop(const GeoLoop *geoloop, const BBox *loopBBox, + const CellBoundary *boundary, + const BBox *boundaryBBox) { + if (!bboxOverlapsBBox(loopBBox, boundaryBBox)) { + return false; + } + LongitudeNormalization loopNormalization; + LongitudeNormalization boundaryNormalization; + bboxNormalization(loopBBox, boundaryBBox, &loopNormalization, + &boundaryNormalization); + + CellBoundary normalBoundary = *boundary; + for (int i = 0; i < boundary->numVerts; i++) { + normalBoundary.verts[i].lng = + normalizeLng(normalBoundary.verts[i].lng, boundaryNormalization); + } + + BBox normalBoundaryBBox = { + .north = boundaryBBox->north, + .south = boundaryBBox->south, + .east = normalizeLng(boundaryBBox->east, boundaryNormalization), + .west = normalizeLng(boundaryBBox->west, boundaryNormalization)}; + + LatLng loop1; + LatLng loop2; + for (int i = 0; i < geoloop->numVerts; i++) { + loop1 = geoloop->verts[i]; + loop1.lng = normalizeLng(loop1.lng, loopNormalization); + loop2 = geoloop->verts[(i + 1) % geoloop->numVerts]; + loop2.lng = normalizeLng(loop2.lng, loopNormalization); + + // Quick check if the line segment overlaps our bbox + if ((loop1.lat >= normalBoundaryBBox.north && + loop2.lat >= normalBoundaryBBox.north) || + (loop1.lat <= normalBoundaryBBox.south && + loop2.lat <= normalBoundaryBBox.south) || + (loop1.lng <= normalBoundaryBBox.west && + loop2.lng <= normalBoundaryBBox.west) || + (loop1.lng >= normalBoundaryBBox.east && + loop2.lng >= normalBoundaryBBox.east)) { + continue; + } + + for (int j = 0; j < normalBoundary.numVerts; j++) { + if (lineIntersectsLine( + &loop1, &loop2, &normalBoundary.verts[j], + &normalBoundary.verts[(j + 1) % normalBoundary.numVerts])) { + return true; + } + } + } + return false; +} + +/** + * Whether two lines intersect. This is a purely Cartesian implementation + * and does not consider anti-meridian wrapping, poles, etc. Based on + * http://www.jeffreythompson.org/collision-detection/line-line.php + * @param a1 Start of line A + * @param a2 End of line A + * @param b1 Start of line B + * @param b2 End of line B + * @return Whether the lines intersect + */ +bool lineIntersectsLine(const LatLng *a1, const LatLng *a2, const LatLng *b1, + const LatLng *b2) { + double denom = ((b2->lng - b1->lng) * (a2->lat - a1->lat) - + (b2->lat - b1->lat) * (a2->lng - a1->lng)); + if (!denom) return false; + + double test; + test = ((b2->lat - b1->lat) * (a1->lng - b1->lng) - + (b2->lng - b1->lng) * (a1->lat - b1->lat)) / + denom; + if (test < 0 || test > 1) return false; + + test = ((a2->lat - a1->lat) * (a1->lng - b1->lng) - + (a2->lng - a1->lng) * (a1->lat - b1->lat)) / + denom; + return (test >= 0 && test <= 1); +}