diff --git a/CHANGELOG.md b/CHANGELOG.md index f078191fde..2fd6506317 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ The public API of this library consists of the functions declared in file - `reverseDirectedEdge` function (#1098) - Bidirectional `gridPathCells` (#1111) - `h3.pc` pkg-config file (#1115) +- Geodesic coverage flag for `polygonToCellsExperimental` that follows great-circle edges when covering very large polygons. This mode is significantly slower than the planar algorithm, so prefer lower resolutions when using it. ### Changed - `cellsToLinkedMultiPolygon` now uses `cellsToMultiPolygon` internally (#1126, #1133, #1136) diff --git a/CMakeLists.txt b/CMakeLists.txt index 92f29201d6..3175d86dd8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -172,6 +172,11 @@ set(LIB_SOURCE_FILES src/h3lib/include/adder.h src/h3lib/include/area.h src/h3lib/include/cellsToMultiPoly.h + src/h3lib/include/polyfillIterator.h + src/h3lib/include/geodesicIterator.h + src/h3lib/include/geodesicPolygonInternal.h + src/h3lib/include/geodesicCellBoundary.h + src/h3lib/include/sphereCapTables.h src/h3lib/lib/h3Assert.c src/h3lib/lib/algos.c src/h3lib/lib/bbox.c @@ -190,7 +195,9 @@ set(LIB_SOURCE_FILES src/h3lib/lib/faceijk.c src/h3lib/lib/baseCells.c src/h3lib/lib/area.c - src/h3lib/lib/cellsToMultiPoly.c) + src/h3lib/lib/cellsToMultiPoly.c + src/h3lib/lib/geodesicIterator.c + src/h3lib/lib/geodesicPolygon.c) set(APP_SOURCE_FILES src/apps/applib/include/kml.h src/apps/applib/include/benchmark.h @@ -222,6 +229,10 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testPolygonToCellsExperimental.c src/apps/testapps/testPolygonToCellsReported.c src/apps/testapps/testPolygonToCellsReportedExperimental.c + src/apps/testapps/testGeodesicPolygonInternal.c + src/apps/testapps/testGeodesicPolygonToCellsExperimental.c + src/apps/testapps/testGeodesicCoverage.c + src/apps/testapps/testGeodesicGridDiskRoundtrip.c src/apps/testapps/testPentagonIndexes.c src/apps/testapps/testGridDisk.c src/apps/testapps/testGridDiskInternal.c @@ -243,6 +254,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testGetIcosahedronFaces.c src/apps/testapps/testLatLng.c src/apps/testapps/testLatLngInternal.c + src/apps/testapps/testVec3d.c src/apps/testapps/testGridRing.c src/apps/testapps/testGridRingInternal.c src/apps/testapps/testGridRingUnsafe.c @@ -284,6 +296,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testDescribeH3Error.c src/apps/testapps/testGosperIter.c src/apps/testapps/testGeoLoopArea.c + src/apps/testapps/testSphereCap.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c @@ -309,6 +322,7 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerLocalIj.c src/apps/fuzzers/fuzzerPolygonToCells.c src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c + src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c src/apps/fuzzers/fuzzerPolygonToCellsNoHoles.c src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c src/apps/fuzzers/fuzzerCellToChildPos.c @@ -316,6 +330,7 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerInternalCoordIjk.c src/apps/benchmarks/benchmarkPolygonToCells.c src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c + src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c src/apps/benchmarks/benchmarkPolygon.c src/apps/benchmarks/benchmarkCellsToPolyAlgos.c src/apps/benchmarks/benchmarkCellToChildren.c @@ -644,6 +659,8 @@ if(BUILD_FUZZERS) add_h3_fuzzer(fuzzerPolygonToCells src/apps/fuzzers/fuzzerPolygonToCells.c) add_h3_fuzzer(fuzzerPolygonToCellsExperimental src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c) + add_h3_fuzzer(fuzzerPolygonToCellsExperimentalGeodesic + src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c) add_h3_fuzzer(fuzzerPolygonToCellsNoHoles src/apps/fuzzers/fuzzerPolygonToCellsNoHoles.c) add_h3_fuzzer(fuzzerPolygonToCellsExperimentalNoHoles @@ -690,6 +707,8 @@ if(BUILD_BENCHMARKS) src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c) add_h3_benchmark(benchmarkArea src/apps/benchmarks/benchmarkArea.c) + add_h3_benchmark(benchmarkPolygonToCellsExperimentalGeodesic + src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_benchmark(benchmarkPolygon src/apps/benchmarks/benchmarkPolygon.c) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index f48c8f9315..32d54f4515 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -243,6 +243,7 @@ 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(testVec3 src/apps/testapps/testVec3.c) +add_h3_test(testVec3d src/apps/testapps/testVec3d.c) add_h3_test(testCellToLocalIj src/apps/testapps/testCellToLocalIj.c) add_h3_test(testCellToLocalIjInternal src/apps/testapps/testCellToLocalIjInternal.c) @@ -270,6 +271,12 @@ add_h3_test_with_arg(testH3NeighborRotations add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 2) +add_h3_test(testGeodesicPolygonToCellsExperimental src/apps/testapps/testGeodesicPolygonToCellsExperimental.c) +add_h3_test(testGeodesicPolygonInternal src/apps/testapps/testGeodesicPolygonInternal.c) +add_h3_test(testGeodesicCoverage src/apps/testapps/testGeodesicCoverage.c) +add_h3_test(testGeodesicGridDiskRoundtrip src/apps/testapps/testGeodesicGridDiskRoundtrip.c) +add_h3_test(testSphereCap src/apps/testapps/testSphereCap.c) + # The "Exhaustive" part of the test name is used by the test-fast to exclude # these files. test-fast exists so that Travis CI can run Valgrind on tests # without taking a very long time. diff --git a/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c b/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c index 9b69dfeb5c..142937fbe3 100644 --- a/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c +++ b/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c @@ -13,6 +13,9 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +/** @file benchmarkPolygonToCellsExperimental.c + * @brief Benchmarks the experimental polygon-to-cells traversal. + */ #include "algos.h" #include "benchmark.h" #include "h3api.h" diff --git a/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c b/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c new file mode 100644 index 0000000000..2bbf89fa2b --- /dev/null +++ b/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c @@ -0,0 +1,208 @@ +/* + * Copyright 2026 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 benchmarkPolygonToCellsExperimentalGeodesic.c + * @brief Benchmarks the experimental polygon-to-cells traversal (geodesic). + */ +#include +#include + +#include "algos.h" +#include "benchmark.h" +#include "h3api.h" +#include "polyfill.h" +#include "polygon.h" + +// Colorado state +LatLng stateVerts[] = {{0.645778804554910, -1.903190792178713}, + {0.645682811446050, -1.780975856637062}, + {0.715595465293187, -1.781167842854781}, + {0.715578012000667, -1.903262350678044}, + {0.645778804554910, -1.903190792178713}}; +GeoLoop stateGeoLoop; +GeoPolygon stateGeoPolygon; + +// Large ellipse approx 1800x1100km +LatLng largeEllipseVerts[] = {{0.698131700797732, -1.658062789394613}, + {0.697959500367918, -1.647103776346370}, + {0.697443578674904, -1.636188013515931}, + {0.696585971825909, -1.625358580432251}, + {0.695390064403362, -1.614658215920214}, + {0.693860576107503, -1.604129149430017}, + {0.692003543129878, -1.593812934376823}, + {0.689826294331253, -1.583750284148419}, + {0.687337422317945, -1.573980911428067}, + {0.684546749530733, -1.564543371466692}, + {0.681465289480172, -1.555474909922921}, + {0.678105203281295, -1.546811315871481}, + {0.674479751659253, -1.538586780560073}, + {0.670603242615285, -1.530833762472139}, + {0.666490974959581, -1.523582859228054}, + {0.662159177933861, -1.516862686830300}, + {0.657624947161976, -1.510699766729177}, + {0.652906177181289, -1.505118421154754}, + {0.648021490821112, -1.500140677128138}, + {0.642990165706910, -1.495786179530887}, + {0.637832058180313, -1.492072113575638}, + {0.632567524935215, -1.489013136983919}, + {0.627217342679196, -1.486621322138826}, + {0.621802626137356, -1.484906108440836}, + {0.616344744722137, -1.483874265054807}, + {0.610865238198015, -1.483529864195180}, + {0.605385731673894, -1.483874265054807}, + {0.599927850258674, -1.484906108440836}, + {0.594513133716834, -1.486621322138826}, + {0.589162951460816, -1.489013136983919}, + {0.583898418215717, -1.492072113575638}, + {0.578740310689120, -1.495786179530887}, + {0.573708985574918, -1.500140677128138}, + {0.568824299214742, -1.505118421154754}, + {0.564105529234055, -1.510699766729177}, + {0.559571298462169, -1.516862686830300}, + {0.555239501436449, -1.523582859228054}, + {0.551127233780745, -1.530833762472139}, + {0.547250724736778, -1.538586780560073}, + {0.543625273114736, -1.546811315871481}, + {0.540265186915859, -1.555474909922921}, + {0.537183726865297, -1.564543371466692}, + {0.534393054078086, -1.573980911428067}, + {0.531904182064778, -1.583750284148419}, + {0.529726933266152, -1.593812934376823}, + {0.527869900288528, -1.604129149430017}, + {0.526340411992668, -1.614658215920214}, + {0.525144504570122, -1.625358580432251}, + {0.524286897721127, -1.636188013515931}, + {0.523770976028112, -1.647103776346370}, + {0.523598775598299, -1.658062789394613}, + {0.523770976028112, -1.669021802442856}, + {0.524286897721127, -1.679937565273295}, + {0.525144504570122, -1.690766998356975}, + {0.526340411992668, -1.701467362869012}, + {0.527869900288528, -1.711996429359209}, + {0.529726933266152, -1.722312644412403}, + {0.531904182064778, -1.732375294640808}, + {0.534393054078086, -1.742144667361160}, + {0.537183726865297, -1.751582207322534}, + {0.540265186915859, -1.760650668866305}, + {0.543625273114736, -1.769314262917745}, + {0.547250724736778, -1.777538798229153}, + {0.551127233780745, -1.785291816317088}, + {0.555239501436449, -1.792542719561173}, + {0.559571298462169, -1.799262891958926}, + {0.564105529234055, -1.805425812060049}, + {0.568824299214742, -1.811007157634472}, + {0.573708985574918, -1.815984901661088}, + {0.578740310689120, -1.820339399258339}, + {0.583898418215717, -1.824053465213588}, + {0.589162951460816, -1.827112441805307}, + {0.594513133716834, -1.829504256650400}, + {0.599927850258674, -1.831219470348391}, + {0.605385731673894, -1.832251313734419}, + {0.610865238198015, -1.832595714594046}, + {0.616344744722137, -1.832251313734419}, + {0.621802626137356, -1.831219470348391}, + {0.627217342679196, -1.829504256650400}, + {0.632567524935215, -1.827112441805307}, + {0.637832058180313, -1.824053465213588}, + {0.642990165706910, -1.820339399258339}, + {0.648021490821112, -1.815984901661088}, + {0.652906177181289, -1.811007157634472}, + {0.657624947161976, -1.805425812060049}, + {0.662159177933861, -1.799262891958926}, + {0.666490974959581, -1.792542719561173}, + {0.670603242615285, -1.785291816317088}, + {0.674479751659253, -1.777538798229153}, + {0.678105203281295, -1.769314262917745}, + {0.681465289480172, -1.760650668866305}, + {0.684546749530733, -1.751582207322534}, + {0.687337422317945, -1.742144667361160}, + {0.689826294331253, -1.732375294640808}, + {0.692003543129878, -1.722312644412403}, + {0.693860576107503, -1.711996429359210}, + {0.695390064403362, -1.701467362869012}, + {0.696585971825909, -1.690766998356975}, + {0.697443578674904, -1.679937565273295}, + {0.697959500367918, -1.669021802442856}}; +GeoLoop largeEllipseGeoLoop; +GeoPolygon largeEllipseGeoPolygon; + +// London to NY flight +LatLng londonNyVerts[] = {{0.8989737191417, -0.0022305307840}, + {0.7105724077059, -1.2916483662309}, + {0.8989737191417, -0.0022305307840}}; +GeoLoop londonNyGeoLoop; +GeoPolygon londonNyGeoPolygon; + +void polygonToCellsBenchmark(const GeoPolygon *polygon, int resolution, + bool geodesic) { + uint32_t flags = CONTAINMENT_OVERLAPPING; + if (geodesic) { + FLAG_SET_GEODESIC(flags); + } + + int64_t numHexagons = 0; + H3_EXPORT(maxPolygonToCellsSizeExperimental) + (polygon, resolution, flags, &numHexagons); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + H3_EXPORT(polygonToCellsExperimental) + (polygon, resolution, flags, numHexagons, hexagons); + free(hexagons); +} + +BEGIN_BENCHMARKS(); + +stateGeoLoop.numVerts = 5; +stateGeoLoop.verts = stateVerts; +stateGeoPolygon.geoloop = stateGeoLoop; + +largeEllipseGeoLoop.numVerts = 100; +largeEllipseGeoLoop.verts = largeEllipseVerts; +largeEllipseGeoPolygon.geoloop = largeEllipseGeoLoop; + +londonNyGeoLoop.numVerts = 3; +londonNyGeoLoop.verts = londonNyVerts; +londonNyGeoPolygon.geoloop = londonNyGeoLoop; + +const int stateResolution = 5; +const int largeEllipseResolution = 4; +const int londonNyResolution = 3; + +BENCHMARK(polygonToCellsState_PlanarOverlapping, 30, { + polygonToCellsBenchmark(&stateGeoPolygon, stateResolution, false); +}); + +BENCHMARK(polygonToCellsState_GeodesicOverlapping, 30, { + polygonToCellsBenchmark(&stateGeoPolygon, stateResolution, true); +}); + +BENCHMARK(polygonToCellsLargeEllipse_PlanarOverlapping, 50, { + polygonToCellsBenchmark(&largeEllipseGeoPolygon, largeEllipseResolution, + false); +}); + +BENCHMARK(polygonToCellsLargeEllipse_GeodesicOverlapping, 50, { + polygonToCellsBenchmark(&largeEllipseGeoPolygon, largeEllipseResolution, + true); +}); + +BENCHMARK(polygonToCellsLondonNY_PlanarOverlapping, 100, { + polygonToCellsBenchmark(&londonNyGeoPolygon, londonNyResolution, false); +}); + +BENCHMARK(polygonToCellsLondonNY_GeodesicOverlapping, 100, { + polygonToCellsBenchmark(&londonNyGeoPolygon, londonNyResolution, true); +}); + +END_BENCHMARKS(); diff --git a/src/apps/fuzzers/README.md b/src/apps/fuzzers/README.md index a644d18533..6a13339a60 100644 --- a/src/apps/fuzzers/README.md +++ b/src/apps/fuzzers/README.md @@ -66,7 +66,7 @@ The public API of H3 is covered in the following fuzzers: | localIjToCell | [fuzzerLocalIj](./fuzzerLocalIj.c) | originToDirectedEdges | [fuzzerDirectedEdge](./fuzzerDirectedEdge.c) | polygonToCells | [fuzzerPoylgonToCells](./fuzzerPolygonToCells.c) -| polygonToCellsExperimental | [fuzzerPoylgonToCellsExperimental](./fuzzerPolygonToCellsExperimental.c) [fuzzerPoylgonToCellsExperimentalNoHoles](./fuzzerPolygonToCellsExperimentalNoHoles.c) +| polygonToCellsExperimental | [fuzzerPoylgonToCellsExperimental](./fuzzerPolygonToCellsExperimental.c) [fuzzerPoylgonToCellsExperimentalGeodesic](./fuzzerPolygonToCellsExperimentalGeodesic.c) [fuzzerPoylgonToCellsExperimentalNoHoles](./fuzzerPolygonToCellsExperimentalNoHoles.c) | radsToDegs | Trivial | stringToH3 | [fuzzerIndexIO](./fuzzerIndexIO.c) | uncompactCells | [fuzzerCompact](./fuzzerCompact.c) diff --git a/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c b/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c index aaa7e4673e..a1302940f2 100644 --- a/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c +++ b/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c @@ -13,8 +13,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -/** @file - * @brief Fuzzer program for polygonToCells2 and related functions +/** @file fuzzerPolygonToCellsExperimental.c + * @brief Fuzzes the experimental polygon-to-cells implementation. */ #include "aflHarness.h" @@ -95,6 +95,12 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { for (uint32_t flags = 0; flags < CONTAINMENT_INVALID; flags++) { geoPolygon.numHoles = originalNumHoles; run(&geoPolygon, flags, res); + + // If polygon has no holes - it has already been tested + // Don't test it again + if (originalNumHoles == 0) { + continue; + } geoPolygon.numHoles = 0; run(&geoPolygon, flags, res); } diff --git a/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c new file mode 100644 index 0000000000..388f01f8f7 --- /dev/null +++ b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c @@ -0,0 +1,121 @@ +/* + * Copyright 2026 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 fuzzerPolygonToCellsExperimentalGeodesic.c + * @brief Fuzzes the experimental polygon-to-cells implementation in geodesic + * mode. + */ + +#include "aflHarness.h" +#include "h3api.h" +#include "polyfill.h" +#include "polygon.h" +#include "utility.h" + +typedef struct { + int res; + int numHoles; + // repeating: num verts, verts + // We add a large fixed buffer so our test case generator for AFL + // knows how large to make the file. + uint8_t buffer[1024]; +} inputArgs; + +const int MAX_GEODESIC_RES = 4; +const int MAX_SZ = 100000; +const int MAX_HOLES = 100; + +static int populateGeoLoop(GeoLoop *g, const uint8_t *data, size_t *offset, + size_t size) { + if (size < *offset + sizeof(int)) { + return 1; + } + int numVerts = *(const int *)(data + *offset); + *offset = *offset + sizeof(int); + g->numVerts = numVerts; + if (size < *offset + sizeof(LatLng) * numVerts) { + return 1; + } + g->verts = (LatLng *)(data + *offset); + *offset = *offset + sizeof(LatLng) * numVerts; + return 0; +} + +static void runGeodesic(GeoPolygon *geoPolygon, uint32_t flags, int res) { + uint32_t geodesicFlags = flags; + FLAG_SET_GEODESIC(geodesicFlags); + + int64_t sz; + H3Error err = H3_EXPORT(maxPolygonToCellsSizeExperimental)( + geoPolygon, res, geodesicFlags, &sz); + if (!err && sz < MAX_SZ) { + H3Index *out = calloc(sz, sizeof(H3Index)); + H3_EXPORT(polygonToCellsExperimental) + (geoPolygon, res, geodesicFlags, sz, out); + free(out); + } +} + +int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { + if (size < sizeof(inputArgs)) { + return 0; + } + const inputArgs *args = (const inputArgs *)data; + int res = args->res % (MAX_GEODESIC_RES + 1); + + GeoPolygon geoPolygon; + int originalNumHoles = args->numHoles % MAX_HOLES; + geoPolygon.numHoles = originalNumHoles; + if (geoPolygon.numHoles < 0) { + return 0; + } + geoPolygon.holes = calloc(geoPolygon.numHoles, sizeof(GeoLoop)); + size_t offset = sizeof(inputArgs) - sizeof(args->buffer); + if (populateGeoLoop(&geoPolygon.geoloop, data, &offset, size)) { + free(geoPolygon.holes); + return 0; + } + for (int i = 0; i < geoPolygon.numHoles; i++) { + if (populateGeoLoop(&geoPolygon.holes[i], data, &offset, size)) { + free(geoPolygon.holes); + return 0; + } + } + + for (uint32_t flags = 0; flags < CONTAINMENT_INVALID; flags++) { + bool canRunGeodesic = + flags == CONTAINMENT_OVERLAPPING || flags == CONTAINMENT_FULL; + canRunGeodesic &= geoPolygon.geoloop.numVerts <= 256; + + if (!canRunGeodesic) { + continue; + } + + geoPolygon.numHoles = originalNumHoles; + runGeodesic(&geoPolygon, flags, res); + + if (originalNumHoles == 0) { + continue; + } + + geoPolygon.numHoles = 0; + runGeodesic(&geoPolygon, flags, res); + } + free(geoPolygon.holes); + + return 0; +} + +AFL_HARNESS_MAIN(sizeof(inputArgs)); diff --git a/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c index d6d7c7e6c4..fcd6faad55 100644 --- a/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c +++ b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c @@ -13,9 +13,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -/** @file - * @brief Fuzzer program for polygonToCellsExperimental and related functions, - * without holes +/** @file fuzzerPolygonToCellsExperimentalNoHoles.c + * @brief Fuzzes the experimental polygon-to-cells path without holes. */ #include "aflHarness.h" diff --git a/src/apps/testapps/testBBoxInternal.c b/src/apps/testapps/testBBoxInternal.c index 15068626cc..8072c860ba 100644 --- a/src/apps/testapps/testBBoxInternal.c +++ b/src/apps/testapps/testBBoxInternal.c @@ -14,6 +14,10 @@ * limitations under the License. */ +/** @file testBBoxInternal.c + * @brief Tests the internal bounding box helpers used by polyfill. + */ + #include #include @@ -23,6 +27,7 @@ #include "polygon.h" #include "test.h" #include "utility.h" +#include "vec3d.h" void assertBBoxFromGeoLoop(const GeoLoop *geoloop, const BBox *expected, const LatLng *inside, const LatLng *outside) { @@ -46,6 +51,13 @@ void assertBBox(const BBox *bbox, const BBox *expected) { t_assert(geoAlmostEqual(&actualSW, &expectedSW), "SW corner matches"); } +bool aabbContainsVec3dWithTolerance(const AABB *box, const Vec3d *v) { + const double tol = 1e-12; + return v->x >= box->min.x - tol && v->x <= box->max.x + tol && + v->y >= box->min.y - tol && v->y <= box->max.y + tol && + v->z >= box->min.z - tol && v->z <= box->max.z + tol; +} + SUITE(BBox) { TEST(posLatPosLng) { LatLng verts[] = {{0.8, 0.3}, {0.7, 0.6}, {1.1, 0.7}, {1.0, 0.2}}; @@ -441,4 +453,86 @@ SUITE(BBox) { scaleBBox(&bbox, 2); assertBBox(&bbox, &expected); } + + TEST(aabbEmptyInvertedState) { + AABB box; + aabbEmptyInverted(&box); + t_assert(box.min.x == 1.0 && box.min.y == 1.0 && box.min.z == 1.0, + "Inverted box initializes min values to 1"); + t_assert(box.max.x == -1.0 && box.max.y == -1.0 && box.max.z == -1.0, + "Inverted box initializes max values to -1"); + } + + TEST(aabbVec3Updates) { + AABB box; + aabbEmptyInverted(&box); + + Vec3d first = {0.5, -0.25, 0.75}; + Vec3d second = {-0.25, 0.5, -0.5}; + + aabbUpdateWithVec3d(&box, &first); + t_assert(box.min.x == 0.5 && box.max.x == 0.5, + "Single point updates x"); + t_assert(box.min.y == -0.25 && box.max.y == -0.25, + "Single point updates y"); + t_assert(box.min.z == 0.75 && box.max.z == 0.75, + "Single point updates z"); + + aabbUpdateWithVec3d(&box, &second); + t_assert(box.min.x == -0.25 && box.max.x == 0.5, + "Second point expands x range"); + t_assert(box.min.y == -0.25 && box.max.y == 0.5, + "Second point expands y range"); + t_assert(box.min.z == -0.5 && box.max.z == 0.75, + "Second point expands z range"); + } + + TEST(aabbArcExtrema) { + const LatLng ll1 = {-0.15745782909055106, -1.4432851587832365}; + const LatLng ll2 = {0.36335746078285036, -1.680739052834902}; + Vec3d v1 = latLngToVec3(ll1); + Vec3d v2 = latLngToVec3(ll2); + Vec3d normal = vec3Cross(v1, v2); + + AABB box; + aabbEmptyInverted(&box); + aabbUpdateWithVec3d(&box, &v1); + aabbUpdateWithVec3d(&box, &v2); + + double minYBefore = box.min.y; + aabbUpdateWithArcExtrema(&box, &v1, &v2, &normal); + t_assert(box.min.y < minYBefore, + "Arc extrema expands bounding box for intermediate point"); + } + + TEST(aabbArcExtremaBoundsTiltedArc) { + Vec3d normal = {1.0, 1.0, 1.0}; + vec3Normalize(&normal); + + Vec3d basis1 = {1.0, -1.0, 0.0}; + vec3Normalize(&basis1); + Vec3d basis2 = vec3Cross(normal, basis1); + vec3Normalize(&basis2); + + const double start = -1.25; + const double end = 1.25; + Vec3d v1 = vec3LinComb(cos(start), basis1, sin(start), basis2); + Vec3d v2 = vec3LinComb(cos(end), basis1, sin(end), basis2); + Vec3d edgeCross = vec3Cross(v1, v2); + + AABB box; + aabbEmptyInverted(&box); + aabbUpdateWithVec3d(&box, &v1); + aabbUpdateWithVec3d(&box, &v2); + aabbUpdateWithArcExtrema(&box, &v1, &v2, &edgeCross); + + for (int i = 0; i <= 128; i++) { + double t = (double)i / 128.0; + double angle = start + (end - start) * t; + Vec3d sample = vec3LinComb(cos(angle), basis1, sin(angle), basis2); + + t_assert(aabbContainsVec3dWithTolerance(&box, &sample), + "Arc extrema AABB contains sampled great-circle point"); + } + } } diff --git a/src/apps/testapps/testGeodesicCoverage.c b/src/apps/testapps/testGeodesicCoverage.c new file mode 100644 index 0000000000..3f7f6a5c59 --- /dev/null +++ b/src/apps/testapps/testGeodesicCoverage.c @@ -0,0 +1,109 @@ +/* + * Copyright 2026 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 testGeodesicCoverage.c + * @brief Additional tests to improve coverage of geodesic functions + */ + +#include + +#include "geodesicPolygonInternal.h" +#include "latLng.h" +#include "polygon.h" +#include "test.h" + +#define DEG_TO_RAD (M_PI / 180.0) + +SUITE(GeodesicCoverage) { + TEST(parallelEdges) { + // Create a polygon with nearly parallel edges to trigger + // edge intersection swap code paths + LatLng verts[] = {{.lat = 0.0, .lng = 0.0}, + {.lat = 0.0, .lng = 1.0 * DEG_TO_RAD}, + {.lat = 0.01 * DEG_TO_RAD, .lng = 1.0 * DEG_TO_RAD}, + {.lat = 0.01 * DEG_TO_RAD, .lng = 0.0}}; + GeoLoop loop = {.numVerts = 4, .verts = verts}; + GeoPolygon polygon = {.geoloop = loop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + // Test boundary intersection with parallel edges + GeodesicCellBoundary boundary = {.numVerts = 2}; + LatLng edgePts[] = { + {.lat = 0.005 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 0.005 * DEG_TO_RAD, .lng = 1.0 * DEG_TO_RAD}}; + for (int i = 0; i < 2; i++) { + boundary.verts[i] = latLngToVec3(edgePts[i]); + } + SphereCap cap = {.cosRadius = cos(2.0 * DEG_TO_RAD)}; + cap.center = latLngToVec3(edgePts[0]); + t_assert(geodesicPolygonBoundaryIntersects(poly, &boundary, &cap), + "boundary crossing near-parallel edges intersects"); + + geodesicPolygonDestroy(poly); + } + + TEST(reversedEdgeProjections) { + // Create edges where projections need to be swapped + // to exercise the projection-order swap path. + LatLng verts[] = {{.lat = 10.0 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 0.0, .lng = 0.0}, + {.lat = 0.0, .lng = 10.0 * DEG_TO_RAD}}; + GeoLoop loop = {.numVerts = 3, .verts = verts}; + GeoPolygon polygon = {.geoloop = loop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + // Create a boundary that might have reversed projections + GeodesicCellBoundary boundary = {.numVerts = 3}; + LatLng boundaryPts[] = { + {.lat = 5.0 * DEG_TO_RAD, .lng = 5.0 * DEG_TO_RAD}, + {.lat = 0.0, .lng = 5.0 * DEG_TO_RAD}, + {.lat = 2.5 * DEG_TO_RAD, .lng = 2.5 * DEG_TO_RAD}}; + for (int i = 0; i < 3; i++) { + boundary.verts[i] = latLngToVec3(boundaryPts[i]); + } + SphereCap cap = {.cosRadius = -0.5}; + cap.center = latLngToVec3(boundaryPts[0]); + t_assert(geodesicPolygonBoundaryIntersects(poly, &boundary, &cap), + "reversed projection boundary intersects"); + + geodesicPolygonDestroy(poly); + } + + TEST(distantAABBSphere) { + // Create a small polygon that doesn't intersect with origin sphere + // to exercise the fast-reject sphere/AABB overlap path. + LatLng verts[] = {{.lat = 89.0 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 89.0 * DEG_TO_RAD, .lng = 1.0 * DEG_TO_RAD}, + {.lat = 89.5 * DEG_TO_RAD, .lng = 0.5 * DEG_TO_RAD}}; + GeoLoop loop = {.numVerts = 3, .verts = verts}; + GeoPolygon polygon = {.geoloop = loop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + // Test with a distant point + LatLng farPt = {.lat = -89.0 * DEG_TO_RAD, .lng = 180.0 * DEG_TO_RAD}; + Vec3d farVec = latLngToVec3(farPt); + t_assert(!geodesicPolygonContainsPoint(poly, &farVec), + "south-pole distant point is outside north-pole polygon"); + + geodesicPolygonDestroy(poly); + } +} diff --git a/src/apps/testapps/testGeodesicGridDiskRoundtrip.c b/src/apps/testapps/testGeodesicGridDiskRoundtrip.c new file mode 100644 index 0000000000..3feaa07026 --- /dev/null +++ b/src/apps/testapps/testGeodesicGridDiskRoundtrip.c @@ -0,0 +1,517 @@ +/* + * Copyright 2026 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 testGeodesicGridDiskRoundtrip.c + * @brief Tests that geodesic polygonToCellsExperimental with + * CONTAINMENT_OVERLAPPING | FLAG_GEODESIC correctly recovers cells + * from polygons generated by cellsToMultiPolygon. + */ + +#include +#include +#include + +#include "cellsToMultiPoly.h" +#include "h3Index.h" +#include "iterators.h" +#include "test.h" +#include "vec3d.h" + +#define TRACE_VEC_MATCH_TOL 1e-12 + +typedef struct { + int missingCells; + int skippedCases; + int internalErrors; +} RoundtripStats; + +typedef struct { + const char *name; + int res; + int k; +} ExhaustiveConfig; + +static double vec3Distance(const Vec3d *a, const Vec3d *b) { + const double dx = a->x - b->x; + const double dy = a->y - b->y; + const double dz = a->z - b->z; + return sqrt(dx * dx + dy * dy + dz * dz); +} + +static int64_t compactNonZeroCells(H3Index *cells, int64_t count) { + int64_t out = 0; + for (int64_t i = 0; i < count; i++) { + if (cells[i] != H3_NULL) { + cells[out++] = cells[i]; + } + } + return out; +} + +static int traceCellBoundaryPathMismatch(H3Index cell) { + CellBoundary latLngBoundary = {0}; + GeodesicCellBoundary geodesicBoundary = {0}; + + if (H3_EXPORT(cellToBoundary)(cell, &latLngBoundary) != E_SUCCESS || + cellToGeodesicBoundary(cell, &geodesicBoundary) != E_SUCCESS) { + fprintf(stderr, + " unable to trace boundary path mismatch for 0x%llx\n", + (unsigned long long)cell); + return 1; + } + + fprintf(stderr, " cell boundary path comparison for 0x%llx:\n", + (unsigned long long)cell); + for (int i = 0; i < latLngBoundary.numVerts; i++) { + Vec3d viaGeo = latLngToVec3(latLngBoundary.verts[i]); + + const double dist = vec3Distance(&viaGeo, &geodesicBoundary.verts[i]); + fprintf(stderr, + " vert[%d]: dist(_geoToVec3d(cellToBoundary), " + "cellToGeodesicBoundary)=%.3e\n", + i, dist); + } + + return 0; +} + +static int traceMissingCellAgainstPolygon(const GeoPolygon *polygon, + H3Index missingCell) { + GeodesicCellBoundary cellBoundary = {0}; + if (cellToGeodesicBoundary(missingCell, &cellBoundary) != E_SUCCESS) { + fprintf(stderr, + " unable to trace missing cell boundary for 0x%llx\n", + (unsigned long long)missingCell); + return 1; + } + + const GeoLoop *loop = &polygon->geoloop; + Vec3d *polyVerts = calloc(loop->numVerts, sizeof(Vec3d)); + if (polyVerts == NULL) { + fprintf(stderr, " alloc failure while tracing missing cell 0x%llx\n", + (unsigned long long)missingCell); + return 1; + } + + for (int i = 0; i < loop->numVerts; i++) { + polyVerts[i] = latLngToVec3(loop->verts[i]); + } + + fprintf(stderr, " nearest polygon vertex for missing cell 0x%llx:\n", + (unsigned long long)missingCell); + for (int i = 0; i < cellBoundary.numVerts; i++) { + int nearestIndex = -1; + double nearestDist = INFINITY; + for (int j = 0; j < loop->numVerts; j++) { + const double dist = + vec3Distance(&cellBoundary.verts[i], &polyVerts[j]); + if (dist < nearestDist) { + nearestDist = dist; + nearestIndex = j; + } + } + + fprintf(stderr, " cell vert[%d] -> poly vert[%d], dist=%.3e\n", i, + nearestIndex, nearestDist); + } + + fprintf(stderr, + " near-coincident edge traces (_geodesicEdgesCross inputs):\n"); + for (int j = 0; j < cellBoundary.numVerts; j++) { + const int nextJ = (j + 1) % cellBoundary.numVerts; + Vec3d normalA = + vec3Cross(cellBoundary.verts[j], cellBoundary.verts[nextJ]); + + for (int i = 0; i < loop->numVerts; i++) { + const int nextI = (i + 1) % loop->numVerts; + + const double d00 = + vec3Distance(&cellBoundary.verts[j], &polyVerts[i]); + const double d11 = + vec3Distance(&cellBoundary.verts[nextJ], &polyVerts[nextI]); + const double d01 = + vec3Distance(&cellBoundary.verts[j], &polyVerts[nextI]); + const double d10 = + vec3Distance(&cellBoundary.verts[nextJ], &polyVerts[i]); + + if (!((d00 < TRACE_VEC_MATCH_TOL && d11 < TRACE_VEC_MATCH_TOL) || + (d01 < TRACE_VEC_MATCH_TOL && d10 < TRACE_VEC_MATCH_TOL))) { + continue; + } + + Vec3d normalB = vec3Cross(polyVerts[i], polyVerts[nextI]); + + const double b1Side = vec3Dot(normalA, polyVerts[i]); + const double b2Side = vec3Dot(normalA, polyVerts[nextI]); + const double a1Side = vec3Dot(normalB, cellBoundary.verts[j]); + const double a2Side = vec3Dot(normalB, cellBoundary.verts[nextJ]); + const double bProd = b1Side * b2Side; + const double aProd = a1Side * a2Side; + const int rejected = (bProd > 0) || (aProd > 0); + + fprintf(stderr, + " poly[%d] x cell[%d]: " + "d00=%.3e d11=%.3e d01=%.3e d10=%.3e " + "b1Side=%.3e b2Side=%.3e a1Side=%.3e a2Side=%.3e\n", + i, j, d00, d11, d01, d10, b1Side, b2Side, a1Side, a2Side); + fprintf(stderr, " -> %s: bProd=%.3e aProd=%.3e\n", + rejected ? "REJECTED (same side)" : "candidate crossing", + bProd, aProd); + } + } + + free(polyVerts); + return 0; +} + +/* ------------------------------------------------------------------ */ +/* Helper: build the outer polygon for gridDisk(center, k). */ +/* Returns E_SUCCESS and populates *mpoly on success. */ +/* Caller must call destroyGeoMultiPolygon on E_SUCCESS. */ +/* ------------------------------------------------------------------ */ +static H3Error buildDiskPolygon(H3Index center, int k, GeoMultiPolygon *mpoly) { + int64_t diskMax = 0; + H3Error err = H3_EXPORT(maxGridDiskSize)(k, &diskMax); + if (err != E_SUCCESS) return err; + + H3Index *disk = calloc((size_t)diskMax, sizeof(H3Index)); + if (!disk) return E_MEMORY_ALLOC; + + err = H3_EXPORT(gridDisk)(center, k, disk); + if (err != E_SUCCESS) { + free(disk); + return err; + } + + const int64_t diskCount = compactNonZeroCells(disk, diskMax); + err = H3_EXPORT(cellsToMultiPolygon)(disk, diskCount, mpoly); + free(disk); + return err; +} + +/* ------------------------------------------------------------------ */ +/* Helper: run the roundtrip check for one (cell, k) pair. */ +/* Updates aggregate stats instead of asserting internally. */ +/* ------------------------------------------------------------------ */ +static void gridDiskRoundtripCheck(H3Index center, int k, + RoundtripStats *stats) { + int res = H3_GET_RESOLUTION(center); + + /* 1. Build grid disk ------------------------------------------- */ + int64_t diskMax; + if (H3_EXPORT(maxGridDiskSize)(k, &diskMax) != E_SUCCESS) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - maxGridDiskSize failed\n", + (unsigned long long)center, k, res); + stats->internalErrors++; + return; + } + H3Index *disk = calloc(diskMax, sizeof(H3Index)); + if (disk == NULL) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - alloc disk failed\n", + (unsigned long long)center, k, res); + stats->internalErrors++; + return; + } + if (H3_EXPORT(gridDisk)(center, k, disk) != E_SUCCESS) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - gridDisk failed\n", + (unsigned long long)center, k, res); + free(disk); + stats->internalErrors++; + return; + } + const int64_t diskCount = compactNonZeroCells(disk, diskMax); + + /* 2. Cells -> polygon ------------------------------------------ */ + GeoMultiPolygon mpoly = {0}; + if (H3_EXPORT(cellsToMultiPolygon)(disk, diskCount, &mpoly) != E_SUCCESS) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - cellsToMultiPolygon " + "failed\n", + (unsigned long long)center, k, res); + free(disk); + stats->internalErrors++; + return; + } + + if (mpoly.numPolygons != 1 || mpoly.polygons[0].numHoles != 0) { + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + free(disk); + stats->skippedCases++; + return; + } + + /* 3. Geodesic polyfill ----------------------------------------- */ + uint32_t flags = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(flags); + + int64_t outMax = 0; + if (H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &mpoly.polygons[0], res, flags, &outMax) != E_SUCCESS) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - " + "maxPolygonToCellsSizeExperimental failed\n", + (unsigned long long)center, k, res); + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + free(disk); + stats->internalErrors++; + return; + } + + H3Index *out = calloc(outMax, sizeof(H3Index)); + if (out == NULL) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - alloc out failed\n", + (unsigned long long)center, k, res); + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + free(disk); + stats->internalErrors++; + return; + } + if (H3_EXPORT(polygonToCellsExperimental)(&mpoly.polygons[0], res, flags, + outMax, out) != E_SUCCESS) { + fprintf(stderr, + " ERROR: center=0x%llx k=%d res=%d - " + "polygonToCellsExperimental failed\n", + (unsigned long long)center, k, res); + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + free(out); + free(disk); + stats->internalErrors++; + return; + } + + /* 4. Check superset property: every disk cell must be in out ---- */ + for (int64_t i = 0; i < diskCount; i++) { + int found = 0; + for (int64_t j = 0; j < outMax; j++) { + if (out[j] == disk[i]) { + found = 1; + break; + } + } + if (!found) { + stats->missingCells++; + fprintf(stderr, + " FAIL: center=0x%llx k=%d res=%d - " + "missing cell 0x%llx\n", + (unsigned long long)center, k, res, + (unsigned long long)disk[i]); + stats->internalErrors += traceCellBoundaryPathMismatch(disk[i]); + stats->internalErrors += + traceMissingCellAgainstPolygon(&mpoly.polygons[0], disk[i]); + } + } + + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + free(out); + free(disk); +} + +static void runExhaustiveConfig(const ExhaustiveConfig *config, + RoundtripStats *stats) { + for (IterCellsResolution it = iterInitRes(config->res); it.h; + iterStepRes(&it)) { + gridDiskRoundtripCheck(it.h, config->k, stats); + } +} + +/* ------------------------------------------------------------------ */ +/* Test suite */ +/* ------------------------------------------------------------------ */ +SUITE(GeodesicGridDiskRoundtrip) { + TEST(exhaustiveCoverageMatrix) { + const ExhaustiveConfig configs[] = {{"res0k0", 0, 0}, {"res0k1", 0, 1}, + {"res0k2", 0, 2}, {"res1k0", 1, 0}, + {"res1k1", 1, 1}, {"res2k0", 2, 0}}; + // {"res1k2", 1, 2}, + // {"res2k0", 2, 0}, {"res2k1", 2, 1}, + // {"res2k2", 2, 2}, {"res3k0", 3, 0}}; + // You are welcome to test more, but that will take time + RoundtripStats total = {0}; + + for (int i = 0; i < (int)(sizeof(configs) / sizeof(configs[0])); i++) { + RoundtripStats stats = {0}; + runExhaustiveConfig(&configs[i], &stats); + total.missingCells += stats.missingCells; + total.skippedCases += stats.skippedCases; + total.internalErrors += stats.internalErrors; + fprintf(stderr, + " exhaustive %s: missing=%d skipped=%d errors=%d\n", + configs[i].name, stats.missingCells, stats.skippedCases, + stats.internalErrors); + } + + t_assert(total.internalErrors == 0, + "exhaustive configs should not hit helper errors"); + t_assert(total.missingCells == 0, + "all exhaustive res/k configs should round-trip"); + } + + // Stage 1: for each res-0 and res-1 cell, build the k=1 disk polygon + // (the union of the cell and its 6 neighbours) and run FULL geodesic + // polyfill. The center cell has no edges on the outer polygon boundary, so + // it must appear in the FULL result. + TEST(ringPolygonFullContainment) { + const int resolutions[] = {0, 1}; + for (int ri = 0; ri < 2; ri++) { + int res = resolutions[ri]; + for (IterCellsResolution it = iterInitRes(res); it.h; + iterStepRes(&it)) { + H3Index center = it.h; + + GeoMultiPolygon mpoly = {0}; + if (buildDiskPolygon(center, 1, &mpoly) != E_SUCCESS) continue; + if (mpoly.numPolygons != 1 || mpoly.polygons[0].numHoles != 0) { + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + continue; + } + + uint32_t flags = CONTAINMENT_FULL; + FLAG_SET_GEODESIC(flags); + + int64_t outMax = 0; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &mpoly.polygons[0], res, flags, &outMax)); + + H3Index *out = + outMax > 0 ? calloc((size_t)outMax, sizeof(H3Index)) : NULL; + t_assert(outMax == 0 || out != NULL, "alloc succeeded"); + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &mpoly.polygons[0], res, flags, outMax, out)); + + bool found = false; + for (int64_t i = 0; i < outMax && !found; i++) { + if (out[i] == center) found = true; + } + t_assert(found, + "center cell must be in FULL result of k=1 disk " + "polygon"); + + free(out); + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + } + } + } + + // Stage 2: same k=1 disk polygon but with a hole punched through the + // center cell. The hole vertices are midpoints between each center-cell + // boundary vertex and the cell center, so the hole lies entirely inside + // the center cell. + // + // Expected behaviour: + // FULL – center cell excluded (its center is inside the hole, so + // point-in-polygon returns false and FULL skips it) + // OVERLAPPING – center cell included (the hole's first vertex maps to + // the center cell, detected by the hole-vertex fallback) + // – all 6 neighbours included (their centers are inside the + // filled region) + TEST(ringPolygonHoleExcludes) { + for (IterCellsResolution it = iterInitRes(0); it.h; iterStepRes(&it)) { + H3Index center = it.h; + int res = H3_GET_RESOLUTION(center); + + GeoMultiPolygon mpoly = {0}; + if (buildDiskPolygon(center, 1, &mpoly) != E_SUCCESS) continue; + if (mpoly.numPolygons != 1 || mpoly.polygons[0].numHoles != 0) { + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + continue; + } + + // Build the k=1 disk cells list to verify OVERLAPPING later. + int64_t diskMax = 0; + t_assertSuccess(H3_EXPORT(maxGridDiskSize)(1, &diskMax)); + H3Index *disk = calloc((size_t)diskMax, sizeof(H3Index)); + t_assert(disk != NULL, "alloc disk succeeded"); + t_assertSuccess(H3_EXPORT(gridDisk)(center, 1, disk)); + const int64_t diskCount = compactNonZeroCells(disk, diskMax); + + // Build hole: midpoints between boundary vertices and cell center, + // with longitude wrapped to [-pi, pi] to handle transmeridian + // cells. + CellBoundary bnd = {0}; + t_assertSuccess(H3_EXPORT(cellToBoundary)(center, &bnd)); + LatLng cellCenter; + t_assertSuccess(H3_EXPORT(cellToLatLng)(center, &cellCenter)); + + LatLng holeVerts[MAX_CELL_BNDRY_VERTS]; + for (int i = 0; i < bnd.numVerts; i++) { + double dlng = bnd.verts[i].lng - cellCenter.lng; + if (dlng > M_PI) dlng -= 2.0 * M_PI; + if (dlng < -M_PI) dlng += 2.0 * M_PI; + holeVerts[i].lat = (bnd.verts[i].lat + cellCenter.lat) / 2.0; + holeVerts[i].lng = cellCenter.lng + dlng / 2.0; + } + GeoLoop holeLoop = {.numVerts = bnd.numVerts, .verts = holeVerts}; + + // Attach the hole to the outer polygon. + GeoPolygon polyWithHole = mpoly.polygons[0]; + polyWithHole.numHoles = 1; + polyWithHole.holes = &holeLoop; + + uint32_t flagsFull = CONTAINMENT_FULL; + FLAG_SET_GEODESIC(flagsFull); + uint32_t flagsOver = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(flagsOver); + + // FULL: center cell must NOT be in result. + int64_t fullMax = 0; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &polyWithHole, res, flagsFull, &fullMax)); + H3Index *fullOut = + fullMax > 0 ? calloc((size_t)fullMax, sizeof(H3Index)) : NULL; + t_assert(fullMax == 0 || fullOut != NULL, "alloc full succeeded"); + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &polyWithHole, res, flagsFull, fullMax, fullOut)); + + bool inFull = false; + for (int64_t i = 0; i < fullMax && !inFull; i++) { + if (fullOut[i] == center) inFull = true; + } + t_assert(!inFull, + "center cell must NOT appear in FULL result when hole " + "is inside it"); + free(fullOut); + + // OVERLAPPING: every k=1 disk cell (including center) must be in + // result. + int64_t overMax = 0; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &polyWithHole, res, flagsOver, &overMax)); + H3Index *overOut = + overMax > 0 ? calloc((size_t)overMax, sizeof(H3Index)) : NULL; + t_assert(overMax == 0 || overOut != NULL, "alloc over succeeded"); + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)( + &polyWithHole, res, flagsOver, overMax, overOut)); + + for (int64_t i = 0; i < diskCount; i++) { + bool found = false; + for (int64_t j = 0; j < overMax && !found; j++) { + if (overOut[j] == disk[i]) found = true; + } + t_assert(found, + "every k=1 disk cell must be in OVERLAPPING result " + "even with hole inside center cell"); + } + free(overOut); + + free(disk); + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + } + } +} diff --git a/src/apps/testapps/testGeodesicPolygonInternal.c b/src/apps/testapps/testGeodesicPolygonInternal.c new file mode 100644 index 0000000000..e64ae2ecd0 --- /dev/null +++ b/src/apps/testapps/testGeodesicPolygonInternal.c @@ -0,0 +1,557 @@ +/* + * Copyright 2026 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 testGeodesicPolygonInternal.c + * @brief Tests the internal geodesic polygon acceleration structures. + */ + +/* + * Focused tests for the geodesic polygon acceleration structures backing the + * geodesic polygonToCellsExperimental flag. + */ + +#include +#include + +#include "geodesicCellBoundary.h" +#include "geodesicPolygonInternal.h" +#include "h3Index.h" +#include "polygon.h" +#include "test.h" +#include "vec3d.h" + +#define DEG_TO_RAD (M_PI / 180.0) + +static LatLng triangleVerts[] = {{.lat = 0.0, .lng = 0.0}, + {.lat = 0.0, .lng = 2.0 * DEG_TO_RAD}, + {.lat = 2.0 * DEG_TO_RAD, .lng = 0.0}}; +static GeoLoop triangleLoop = {.numVerts = 3, .verts = triangleVerts}; + +static GeoPolygon zeroLoopPolygon = { + .geoloop = {.numVerts = 0, .verts = NULL}, .numHoles = 0, .holes = NULL}; + +SUITE(GeodesicPolygonInternal) { + TEST(createAndDestroy) { + GeoPolygon polygon = { + .geoloop = triangleLoop, .numHoles = 0, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + t_assert(poly != NULL, "triangle polygon builds geodesic structure"); + t_assert(poly->numHoles == 0, "no holes copied"); + geodesicPolygonDestroy(poly); + } + + TEST(invalidInputs) { + GeodesicPolygon *poly = NULL; + t_assert(geodesicPolygonCreate(&zeroLoopPolygon, &poly) == E_DOMAIN, + "zero-vertex polygon rejected"); + + GeoLoop invalidHole = {.numVerts = 0, .verts = NULL}; + GeoPolygon polyWithInvalidHole = { + .geoloop = triangleLoop, .numHoles = 1, .holes = &invalidHole}; + t_assert(geodesicPolygonCreate(&polyWithInvalidHole, &poly) == E_DOMAIN, + "hole with zero vertices rejected"); + } + + TEST(containsPoint) { + GeoPolygon polygon = { + .geoloop = triangleLoop, .numHoles = 0, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + LatLng insideLl = {.lat = 0.5 * DEG_TO_RAD, .lng = 0.5 * DEG_TO_RAD}; + LatLng outsideLl = {.lat = 3.0 * DEG_TO_RAD, .lng = 3.0 * DEG_TO_RAD}; + + Vec3d insideVec = latLngToVec3(insideLl); + Vec3d outsideVec = latLngToVec3(outsideLl); + + t_assert(geodesicPolygonContainsPoint(poly, &insideVec), + "point inside polygon detected"); + t_assert(!geodesicPolygonContainsPoint(poly, &outsideVec), + "point outside polygon rejected"); + + geodesicPolygonDestroy(poly); + } + + TEST(capIntersection) { + GeoPolygon polygon = { + .geoloop = triangleLoop, .numHoles = 0, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + LatLng centerLl = {.lat = 0.5 * DEG_TO_RAD, .lng = 0.5 * DEG_TO_RAD}; + H3Index cell; + t_assertSuccess(H3_EXPORT(latLngToCell)(¢erLl, 1, &cell)); + + SphereCap cap; + t_assertSuccess(cellToSphereCap(cell, &cap)); + t_assert(geodesicPolygonCapIntersects(poly, &cap), + "cap overlapping polygon detected"); + + LatLng farLl = {.lat = 30.0 * DEG_TO_RAD, .lng = -50.0 * DEG_TO_RAD}; + t_assertSuccess(H3_EXPORT(latLngToCell)(&farLl, 1, &cell)); + t_assertSuccess(cellToSphereCap(cell, &cap)); + t_assert(!geodesicPolygonCapIntersects(poly, &cap), + "distant caps rejected"); + + geodesicPolygonDestroy(poly); + } + + TEST(boundaryIntersection) { + GeoPolygon polygon = { + .geoloop = triangleLoop, .numHoles = 0, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + GeodesicCellBoundary boundary = {.numVerts = triangleLoop.numVerts}; + for (int i = 0; i < triangleLoop.numVerts; i++) { + boundary.verts[i] = latLngToVec3(triangleLoop.verts[i]); + } + + SphereCap permissiveCap = {.cosRadius = -1.0}; + permissiveCap.center = latLngToVec3(triangleLoop.verts[0]); + t_assert( + geodesicPolygonBoundaryIntersects(poly, &boundary, &permissiveCap), + "coincident boundaries reported as intersecting"); + + GeodesicCellBoundary farBoundary = {.numVerts = 4}; + LatLng squareLl[] = { + {.lat = 10.0 * DEG_TO_RAD, .lng = 10.0 * DEG_TO_RAD}, + {.lat = 10.0 * DEG_TO_RAD, .lng = 15.0 * DEG_TO_RAD}, + {.lat = 15.0 * DEG_TO_RAD, .lng = 15.0 * DEG_TO_RAD}, + {.lat = 15.0 * DEG_TO_RAD, .lng = 10.0 * DEG_TO_RAD}}; + for (int i = 0; i < farBoundary.numVerts; i++) { + farBoundary.verts[i] = latLngToVec3(squareLl[i]); + } + SphereCap farCap = {.cosRadius = cos(2.0 * DEG_TO_RAD)}; + farCap.center = latLngToVec3(squareLl[0]); + t_assert( + !geodesicPolygonBoundaryIntersects(poly, &farBoundary, &farCap), + "far boundary does not intersect"); + + geodesicPolygonDestroy(poly); + } + + TEST(polygonWithHolesNullPointer) { + GeoPolygon polyWithNullHoles = { + .geoloop = triangleLoop, .numHoles = 1, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assert(geodesicPolygonCreate(&polyWithNullHoles, &poly) == E_DOMAIN, + "polygon with holes > 0 but NULL holes pointer rejected"); + } + + TEST(destroyNullPolygon) { + // This should not crash + geodesicPolygonDestroy(NULL); + t_assert(true, "destroying NULL polygon is safe"); + } + + TEST(colinearEdgeIntersection) { + // Test edge case with nearly colinear edges that require swapping + // projections in _edgeIntersectsEdge + GeoPolygon polygon = { + .geoloop = triangleLoop, .numHoles = 0, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + // Test with overlapping boundary that may trigger swap paths + GeodesicCellBoundary boundary = {.numVerts = 3}; + LatLng edgeLl[] = {{.lat = 0.5 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 1.0 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 1.5 * DEG_TO_RAD, .lng = 0.0}}; + for (int i = 0; i < boundary.numVerts; i++) { + boundary.verts[i] = latLngToVec3(edgeLl[i]); + } + + SphereCap cap = {.cosRadius = cos(2.0 * DEG_TO_RAD)}; + cap.center = latLngToVec3(edgeLl[0]); + t_assert(geodesicPolygonBoundaryIntersects(poly, &boundary, &cap), + "colinear overlapping segment intersects polygon boundary"); + + geodesicPolygonDestroy(poly); + } + + TEST(polygonWithAntipodal) { + // Test polygon with edges that may create edge cases in intersection + LatLng antipodeVerts[] = {{.lat = 0.0, .lng = 0.0}, + {.lat = 0.0, .lng = M_PI}, + {.lat = 1.0 * DEG_TO_RAD, .lng = 0.0}}; + GeoLoop antipodeLoop = {.numVerts = 3, .verts = antipodeVerts}; + GeoPolygon polygon = { + .geoloop = antipodeLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + // Opposite-hemisphere point should be rejected quickly. + LatLng oppositePt = {.lat = 0.0, .lng = M_PI}; + Vec3d oppositeVec = latLngToVec3(oppositePt); + t_assert(!geodesicPolygonContainsPoint(poly, &oppositeVec), + "opposite-hemisphere point is outside antipodal polygon"); + + geodesicPolygonDestroy(poly); + } + + TEST(largePolygonContainsCardinalPoints) { + // Test a large polygon that contains cardinal axis points + // to exercise cardinal-axis probing in _geodesicLoopToAABB. + LatLng largeVerts[] = { + {.lat = 60.0 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 60.0 * DEG_TO_RAD, .lng = 90.0 * DEG_TO_RAD}, + {.lat = 60.0 * DEG_TO_RAD, .lng = 180.0 * DEG_TO_RAD}, + {.lat = 60.0 * DEG_TO_RAD, .lng = -90.0 * DEG_TO_RAD}}; + GeoLoop largeLoop = {.numVerts = 4, .verts = largeVerts}; + GeoPolygon largePolygon = { + .geoloop = largeLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&largePolygon, &poly)); + + // Test with north pole which should be inside + LatLng northPole = {.lat = M_PI_2, .lng = 0.0}; + Vec3d northVec = latLngToVec3(northPole); + t_assert(geodesicPolygonContainsPoint(poly, &northVec), + "north pole is inside high-latitude polygon"); + + geodesicPolygonDestroy(poly); + } + + TEST(verySmallPolygon) { + // Test a very small polygon to potentially trigger AABB edge cases + LatLng tinyVerts[] = { + {.lat = 0.00001 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 0.00001 * DEG_TO_RAD, .lng = 0.00001 * DEG_TO_RAD}, + {.lat = 0.0, .lng = 0.00001 * DEG_TO_RAD}}; + GeoLoop tinyLoop = {.numVerts = 3, .verts = tinyVerts}; + GeoPolygon tinyPolygon = { + .geoloop = tinyLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&tinyPolygon, &poly)); + + LatLng farPoint = {.lat = 45.0 * DEG_TO_RAD, .lng = 45.0 * DEG_TO_RAD}; + Vec3d farVec = latLngToVec3(farPoint); + t_assert(!geodesicPolygonContainsPoint(poly, &farVec), + "distant point is outside tiny polygon"); + + // Also test cap intersection with very distant cap + H3Index cell; + t_assertSuccess(H3_EXPORT(latLngToCell)(&farPoint, 5, &cell)); + SphereCap cap; + t_assertSuccess(cellToSphereCap(cell, &cap)); + t_assert(!geodesicPolygonCapIntersects(poly, &cap), + "distant cap does not intersect tiny polygon AABB"); + + geodesicPolygonDestroy(poly); + } + + TEST(greatCircleBoundaryLoopRejected) { + // All vertices lie on the equator (a single great circle), yielding a + // zero-area boundary loop. Reject to avoid ambiguous containment. + LatLng equatorLineVerts[] = {{.lat = 0.0, .lng = -60.0 * DEG_TO_RAD}, + {.lat = 0.0, .lng = -20.0 * DEG_TO_RAD}, + {.lat = 0.0, .lng = 20.0 * DEG_TO_RAD}, + {.lat = 0.0, .lng = 60.0 * DEG_TO_RAD}}; + GeoLoop equatorLineLoop = {.numVerts = 4, .verts = equatorLineVerts}; + GeoPolygon polygon = { + .geoloop = equatorLineLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assert(geodesicPolygonCreate(&polygon, &poly) == E_DOMAIN, + "great-circle boundary loop rejected in geodesic mode"); + geodesicPolygonDestroy(poly); + } + + TEST(nearGreatCircleLoopAccepted) { + // Similar to the great-circle boundary case, but one vertex is a few + // degrees off the equator so the loop encloses non-zero area. + LatLng nearEquatorVerts[] = { + {.lat = 0.0, .lng = -60.0 * DEG_TO_RAD}, + {.lat = 0.0, .lng = -20.0 * DEG_TO_RAD}, + {.lat = 3.0 * DEG_TO_RAD, .lng = 20.0 * DEG_TO_RAD}, + {.lat = 0.0, .lng = 60.0 * DEG_TO_RAD}}; + GeoLoop nearEquatorLoop = {.numVerts = 4, .verts = nearEquatorVerts}; + GeoPolygon polygon = { + .geoloop = nearEquatorLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + t_assert(poly != NULL, "near-great-circle loop remains valid"); + geodesicPolygonDestroy(poly); + } + + TEST(polarCapSmallCircleLoopAccepted) { + // Vertices share a high latitude (small circle, not a great circle). + // This is a valid polygon and should be accepted. + LatLng polarCapVerts[] = { + {.lat = 80.0 * DEG_TO_RAD, .lng = -170.0 * DEG_TO_RAD}, + {.lat = 80.0 * DEG_TO_RAD, .lng = -60.0 * DEG_TO_RAD}, + {.lat = 80.0 * DEG_TO_RAD, .lng = 60.0 * DEG_TO_RAD}, + {.lat = 80.0 * DEG_TO_RAD, .lng = 170.0 * DEG_TO_RAD}}; + GeoLoop polarCapLoop = {.numVerts = 4, .verts = polarCapVerts}; + GeoPolygon polygon = { + .geoloop = polarCapLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + t_assert(poly != NULL, "polar-cap small-circle loop is accepted"); + + LatLng northPole = {.lat = M_PI_2, .lng = 0.0}; + Vec3d northVec = latLngToVec3(northPole); + t_assert(geodesicPolygonContainsPoint(poly, &northVec), + "north pole is inside polar-cap polygon"); + + geodesicPolygonDestroy(poly); + } + + TEST(hemispherePoleFallsBackToVertexCandidate) { + // This loop is crafted so the centroid candidate is rejected, while + // one vertex is a valid containing-hemisphere pole. + LatLng verts[] = {{.lat = 29.512732152343 * DEG_TO_RAD, + .lng = -159.155546799103 * DEG_TO_RAD}, + {.lat = -11.073667682726 * DEG_TO_RAD, + .lng = -154.999100885587 * DEG_TO_RAD}, + {.lat = 73.688809451961 * DEG_TO_RAD, + .lng = -161.574247136150 * DEG_TO_RAD}, + {.lat = -79.396090284140 * DEG_TO_RAD, + .lng = 63.249426974972 * DEG_TO_RAD}}; + GeoLoop loop = {.numVerts = 4, .verts = verts}; + GeoPolygon polygon = {.geoloop = loop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + t_assert(poly->geoloop.hasHemispherePole, "hemisphere pole is stored"); + + Vec3d expectedPole = latLngToVec3(verts[1]); + t_assert(vec3Dot(poly->geoloop.hemispherePole, expectedPole) > 0.9999, + "fallback pole matches the vertex candidate"); + + geodesicPolygonDestroy(poly); + } + + TEST(hemispherePoleFallsBackToOppositeCrossCandidate) { + // This loop is crafted so centroid and vertex candidates are rejected, + // and the accepted pole comes from -cross(v0, v1). + LatLng verts[] = {{.lat = -39.222488834958 * DEG_TO_RAD, + .lng = 17.278671407850 * DEG_TO_RAD}, + {.lat = 77.709767848239 * DEG_TO_RAD, + .lng = -177.944713397859 * DEG_TO_RAD}, + {.lat = 48.221389544616 * DEG_TO_RAD, + .lng = 115.374928293173 * DEG_TO_RAD}, + {.lat = 65.650528740421 * DEG_TO_RAD, + .lng = 86.581228259951 * DEG_TO_RAD}}; + GeoLoop loop = {.numVerts = 4, .verts = verts}; + GeoPolygon polygon = {.geoloop = loop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + t_assert(poly->geoloop.hasHemispherePole, "hemisphere pole is stored"); + + Vec3d v0 = latLngToVec3(verts[0]); + Vec3d v1 = latLngToVec3(verts[1]); + Vec3d expectedPole = vec3Cross(v0, v1); + expectedPole.x = -expectedPole.x; + expectedPole.y = -expectedPole.y; + expectedPole.z = -expectedPole.z; + vec3Normalize(&expectedPole); + + t_assert(vec3Dot(poly->geoloop.hemispherePole, expectedPole) > 0.9999, + "fallback pole matches opposite cross-product candidate"); + + geodesicPolygonDestroy(poly); + } + + TEST(nullArgumentGuards) { + LatLng testLl = {.lat = 1.0 * DEG_TO_RAD, .lng = 1.0 * DEG_TO_RAD}; + Vec3d testVec = latLngToVec3(testLl); + GeodesicPolygon dummy = {0}; + + SphereCap cap = {.center = testVec, .cosRadius = -1.0}; + GeodesicCellBoundary boundary = {.numVerts = 1, .verts = {testVec}}; + + t_assert(!geodesicPolygonCapIntersects(NULL, &cap), + "null polygon rejected for cap test"); + t_assert(!geodesicPolygonCapIntersects(&dummy, NULL), + "null cap rejected"); + + t_assert(!geodesicPolygonBoundaryIntersects(NULL, &boundary, &cap), + "null polygon rejected for boundary test"); + t_assert(!geodesicPolygonBoundaryIntersects(&dummy, NULL, &cap), + "null boundary rejected"); + t_assert(!geodesicPolygonBoundaryIntersects(&dummy, &boundary, NULL), + "null cap rejected for boundary test"); + + t_assert(!geodesicPolygonContainsPoint(NULL, &testVec), + "null polygon rejected for containment test"); + t_assert(!geodesicPolygonContainsPoint(&dummy, NULL), + "null point rejected for containment test"); + } + + TEST(largePolygonGlobalCoverage) { + // Polar cap at 30 deg N with 8 equally spaced vertices. + // This covers approximately 25% of the sphere (the northern cap above + // the geodesic boundary). With approximately uniform sampling of the + // whole globe, expect ~25% of points inside. Assert <=50% to catch + // the case where the algorithm misidentifies the complement region + // (75%) as the interior. + LatLng capVerts[] = { + {.lat = 30.0 * DEG_TO_RAD, .lng = 0.0}, + {.lat = 30.0 * DEG_TO_RAD, .lng = 45.0 * DEG_TO_RAD}, + {.lat = 30.0 * DEG_TO_RAD, .lng = 90.0 * DEG_TO_RAD}, + {.lat = 30.0 * DEG_TO_RAD, .lng = 135.0 * DEG_TO_RAD}, + {.lat = 30.0 * DEG_TO_RAD, .lng = 180.0 * DEG_TO_RAD}, + {.lat = 30.0 * DEG_TO_RAD, .lng = -135.0 * DEG_TO_RAD}, + {.lat = 30.0 * DEG_TO_RAD, .lng = -90.0 * DEG_TO_RAD}, + {.lat = 30.0 * DEG_TO_RAD, .lng = -45.0 * DEG_TO_RAD}}; + GeoLoop capLoop = {.numVerts = 8, .verts = capVerts}; + GeoPolygon capPolygon = { + .geoloop = capLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&capPolygon, &poly)); + + // Sanity: north pole must be inside, south pole outside. + LatLng northPole = {.lat = M_PI_2, .lng = 0.0}; + Vec3d northVec = latLngToVec3(northPole); + t_assert(geodesicPolygonContainsPoint(poly, &northVec), + "north pole is inside cap polygon"); + + LatLng southPole = {.lat = -M_PI_2, .lng = 0.0}; + Vec3d southVec = latLngToVec3(southPole); + t_assert(!geodesicPolygonContainsPoint(poly, &southVec), + "south pole is outside cap polygon"); + + // Sample points approximately uniformly on the sphere. + // Rows are uniform in sin(lat) for area-preserving distribution. + int insideCount = 0; + int totalCount = 0; + const int nRows = 90; + const int nCols = 180; + for (int i = 0; i <= nRows; i++) { + double sinLat = -1.0 + 2.0 * i / nRows; + double lat = asin(sinLat); + for (int j = 0; j < nCols; j++) { + double lng = -M_PI + 2.0 * M_PI * j / nCols; + LatLng ll = {.lat = lat, .lng = lng}; + Vec3d vec = latLngToVec3(ll); + if (geodesicPolygonContainsPoint(poly, &vec)) { + insideCount++; + } + totalCount++; + } + } + + // Polygon covers ~25% of the sphere. Use generous <=50% bound. + double fraction = (double)insideCount / totalCount; + t_assert(fraction <= 0.30, + "cap polygon contains at most 30%% of globe sample points"); + // Ensure some meaningful fraction is inside (not all-reject). + t_assert(fraction > 0.20, + "cap polygon contains at least 20%% of globe sample points"); + + geodesicPolygonDestroy(poly); + } + + TEST(nearHemisphericPolygonGlobalCoverage) { + // Polar cap at 10 deg N with 12 equally spaced vertices. + // Covers approximately 41% of the sphere + LatLng capVerts[12]; + for (int i = 0; i < 12; i++) { + capVerts[i].lat = 10.0 * DEG_TO_RAD; + capVerts[i].lng = (i * 30.0) * DEG_TO_RAD; + } + GeoLoop capLoop = {.numVerts = 12, .verts = capVerts}; + GeoPolygon capPolygon = { + .geoloop = capLoop, .numHoles = 0, .holes = NULL}; + + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&capPolygon, &poly)); + + // Sanity: north pole inside, equatorial point outside. + LatLng northPole = {.lat = M_PI_2, .lng = 0.0}; + Vec3d northVec = latLngToVec3(northPole); + t_assert(geodesicPolygonContainsPoint(poly, &northVec), + "north pole is inside near-hemispheric cap"); + + LatLng equatorPt = {.lat = 0.0, .lng = 0.0}; + Vec3d equatorVec = latLngToVec3(equatorPt); + t_assert(!geodesicPolygonContainsPoint(poly, &equatorVec), + "equator point is outside near-hemispheric cap"); + + // Uniform spherical sampling + int insideCount = 0; + int totalCount = 0; + const int nRows = 90; + const int nCols = 180; + for (int i = 0; i <= nRows; i++) { + double sinLat = -1.0 + 2.0 * i / nRows; + double lat = asin(sinLat); + for (int j = 0; j < nCols; j++) { + double lng = -M_PI + 2.0 * M_PI * j / nCols; + LatLng ll = {.lat = lat, .lng = lng}; + Vec3d vec = latLngToVec3(ll); + if (geodesicPolygonContainsPoint(poly, &vec)) { + insideCount++; + } + totalCount++; + } + } + + double fraction = (double)insideCount / totalCount; + t_assert(fraction <= 0.45, + "near-hemispheric cap contains at most 45%% of globe sample " + "points"); + t_assert(fraction > 0.35, + "near-hemispheric cap contains at least 35%% of globe sample " + "points"); + + geodesicPolygonDestroy(poly); + } + + TEST(capRejectsAabbOutsideUnitSphere) { + GeodesicPolygon *poly = calloc(1, sizeof(GeodesicPolygon)); + t_assert(poly != NULL, "allocated synthetic polygon"); + + poly->aabb.min = (Vec3d){2.0, 2.0, 2.0}; + poly->aabb.max = (Vec3d){3.0, 3.0, 3.0}; + + LatLng centerLl = {.lat = 0.0, .lng = 0.0}; + SphereCap cap = {.cosRadius = -1.0}; + cap.center = latLngToVec3(centerLl); + + t_assert(!geodesicPolygonCapIntersects(poly, &cap), + "AABB outside unit sphere rejected"); + + geodesicPolygonDestroy(poly); + } + + TEST(degenerateBoundarySegmentRejected) { + GeoPolygon polygon = { + .geoloop = triangleLoop, .numHoles = 0, .holes = NULL}; + GeodesicPolygon *poly = NULL; + t_assertSuccess(geodesicPolygonCreate(&polygon, &poly)); + + // One repeated boundary vertex creates a zero-length segment. + GeodesicCellBoundary boundary = {.numVerts = 1}; + boundary.verts[0] = poly->geoloop.edges[0].vert; + + SphereCap cap = {.center = boundary.verts[0], .cosRadius = -1.0}; + t_assert(!geodesicPolygonBoundaryIntersects(poly, &boundary, &cap), + "degenerate boundary segment does not intersect"); + + geodesicPolygonDestroy(poly); + } +} diff --git a/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c new file mode 100644 index 0000000000..252bbfb63e --- /dev/null +++ b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c @@ -0,0 +1,672 @@ +/* + * Copyright 2026 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 testGeodesicPolygonToCellsExperimental.c + * @brief Tests the geodesic polygon-to-cells experimental API surface. + */ + +/* + * Geodesic polygonToCellsExperimental tests. + */ + +#include +#include +#include + +#include "h3Index.h" +#include "iterators.h" +#include "polyfill.h" +#include "polygon.h" +#include "test.h" +#include "utility.h" + +// Fixtures - Large scale polygons for performance testing +// Large polygon covering western United States (approx 1000km x 1000km) +static LatLng sfVerts[] = {{0.698132, -2.094395}, // 40°N, 120°W + {0.698132, -1.919862}, // 40°N, 110°W + {0.610865, -1.919862}, // 35°N, 110°W + {0.610865, -2.094395}, // 35°N, 120°W + {0.628319, -2.181662}, // 36°N, 125°W + {0.663225, -2.181662}, // 38°N, 125°W + {0.680678, -2.138763}, // 39°N, 122.5°W + {0.680678, -2.050796}}; // 39°N, 117.5°W +static GeoLoop sfGeoLoop = {.numVerts = 8, .verts = sfVerts}; +static GeoPolygon sfGeoPolygon; + +// Large hole covering parts of Nevada/Utah (approx 300km x 200km) +static LatLng holeVerts[] = {{0.663225, -2.007129}, // 38°N, 115°W + {0.645772, -2.007129}, // 37°N, 115°W + {0.645772, -1.989676}, // 37°N, 114°W + {0.663225, -1.989676}, // 38°N, 114°W + {0.680678, -1.998403}}; // 39°N, 114.5°W +static GeoLoop holeGeoLoop = {.numVerts = 5, .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 outOfBoundsVert[] = {{-2000, -2000}}; +static GeoLoop outOfBoundsVertGeoLoop = {.numVerts = 1, + .verts = outOfBoundsVert}; +static GeoPolygon outOfBoundsVertGeoPolygon; + +static LatLng invalid2Verts[] = {{NAN, NAN}, {-NAN, -NAN}}; +static GeoLoop invalid2GeoLoop = {.numVerts = 2, .verts = invalid2Verts}; +static GeoPolygon invalid2GeoPolygon; + +static GeoLoop nullGeoLoop = {.numVerts = 0}; +static GeoPolygon nullGeoPolygon; + +static LatLng pointVerts[] = {{0.6595072188743, -2.1371053983433}}; +static GeoLoop pointGeoLoop = {.numVerts = 1, .verts = pointVerts}; +static GeoPolygon pointGeoPolygon; + +static LatLng lineVerts[] = {{0.6595072188743, -2.1371053983433}, + {0.6591482046471, -2.1373141048153}}; +static GeoLoop lineGeoLoop = {.numVerts = 2, .verts = lineVerts}; +static GeoPolygon lineGeoPolygon; + +// Unsupported in geodesic mode: loop cannot fit into a hemisphere. +// Very wide loop that wraps almost all longitudes with low area. +static LatLng wideWrapVerts[] = {{10.0 * M_PI / 180.0, -170.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, -110.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, -50.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, 10.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, 70.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, 130.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, 170.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, 170.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, 130.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, 70.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, 10.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, -50.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, -110.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, -170.0 * M_PI / 180.0}}; +static GeoLoop wideWrapGeoLoop = {.numVerts = ARRAY_SIZE(wideWrapVerts), + .verts = wideWrapVerts}; +static GeoPolygon wideWrapGeoPolygon; + +// Very tall north-south loop crossing multiple hemispheres. +static LatLng tallSliceVerts[] = {{60.0 * M_PI / 180.0, 0.0}, + {0.0, 120.0 * M_PI / 180.0}, + {-60.0 * M_PI / 180.0, 0.0}, + {0.0, -120.0 * M_PI / 180.0}}; +static GeoLoop tallSliceGeoLoop = {.numVerts = ARRAY_SIZE(tallSliceVerts), + .verts = tallSliceVerts}; +static GeoPolygon tallSliceGeoPolygon; + +// Another unsupported loop that cannot fit into one hemisphere. +static LatLng reversedTriangleVerts[] = { + {10.0 * M_PI / 180.0, 0.0}, + {10.0 * M_PI / 180.0, 120.0 * M_PI / 180.0}, + {10.0 * M_PI / 180.0, -120.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, -120.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, 120.0 * M_PI / 180.0}, + {-10.0 * M_PI / 180.0, 0.0}}; +static GeoLoop reversedTriangleGeoLoop = { + .numVerts = ARRAY_SIZE(reversedTriangleVerts), + .verts = reversedTriangleVerts}; +static GeoPolygon reversedTriangleGeoPolygon; + +// Large transmeridian polygon covering Pacific (approx 6000km x 4000km) +static LatLng transVerts[] = {{0.698132, -2.967060}, // 40°N, 170°W + {0.698132, 2.967060}, // 40°N, 170°E + {0.174533, 2.967060}, // 10°N, 170°E + {0.174533, -2.967060}, // 10°N, 170°W + {0.436332, -3.141593}, // 25°N, 180°W + {0.610865, -3.141593}}; // 35°N, 180°W +static GeoLoop transGeoLoop = {.numVerts = ARRAY_SIZE(transVerts), + .verts = transVerts}; +static GeoPolygon transGeoPolygon; + +static GeoPolygon nullHoleGeoPolygon; +static GeoPolygon pointHoleGeoPolygon; +static GeoPolygon lineHoleGeoPolygon; + +static int64_t geodesicFillCount(const GeoPolygon *polygon, int res, + uint32_t mode) { + uint32_t flags = mode; + FLAG_SET_GEODESIC(flags); + + int64_t size = 0; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(polygon, res, + flags, &size)); + + H3Index *cells = size > 0 ? calloc(size, sizeof(H3Index)) : NULL; + t_assert(size == 0 || cells != NULL, "allocated output buffer"); + + t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(polygon, res, flags, + size, cells)); + + int64_t count = countNonNullIndexes(cells, size); + free(cells); + return count; +} + +static H3Error geodesicFill(const GeoPolygon *polygon, int res, uint32_t mode, + int64_t *sizeOut, H3Index **cellsOut) { + uint32_t flags = mode; + FLAG_SET_GEODESIC(flags); + + int64_t size = 0; + H3Error err = H3_EXPORT(maxPolygonToCellsSizeExperimental)(polygon, res, + flags, &size); + if (err != E_SUCCESS) { + return err; + } + + H3Index *cells = size > 0 ? calloc(size, sizeof(H3Index)) : NULL; + if (size > 0) { + t_assert(cells != NULL, "allocated output buffer"); + } + + err = + H3_EXPORT(polygonToCellsExperimental)(polygon, res, flags, size, cells); + if (err != E_SUCCESS) { + free(cells); + return err; + } + + if (sizeOut) { + *sizeOut = size; + } else { + t_assert(size == 0 || cells != NULL, "buffer consistent"); + } + if (cellsOut) { + *cellsOut = cells; + } else if (cells) { + free(cells); + } + return E_SUCCESS; +} + +static int _cmpH3Index(const void *a, const void *b) { + H3Index ha = *(const H3Index *)a; + H3Index hb = *(const H3Index *)b; + if (ha < hb) return -1; + if (ha > hb) return 1; + return 0; +} + +// Fill polygon geodesically, compact non-null results and sort them in-place. +// Returns the count of non-null results, or -1 on error. Caller must free *out. +static int64_t _geodesicFillSorted(const GeoPolygon *polygon, int res, + uint32_t mode, H3Index **out) { + uint32_t flags = mode; + FLAG_SET_GEODESIC(flags); + + int64_t sz = 0; + if (H3_EXPORT(maxPolygonToCellsSizeExperimental)(polygon, res, flags, + &sz) != E_SUCCESS) { + return -1; + } + + H3Index *cells = sz > 0 ? calloc((size_t)sz, sizeof(H3Index)) : NULL; + if (sz > 0 && !cells) return -1; + + if (H3_EXPORT(polygonToCellsExperimental)(polygon, res, flags, sz, cells) != + E_SUCCESS) { + free(cells); + return -1; + } + + int64_t n = 0; + for (int64_t i = 0; i < sz; i++) { + if (cells[i] != H3_NULL) cells[n++] = cells[i]; + } + qsort(cells, (size_t)n, sizeof(H3Index), _cmpH3Index); + + *out = cells; + return n; +} + +// Whether `target` appears in a cell array. +static bool _containsCell(const H3Index *cells, int64_t n, H3Index target) { + for (int64_t i = 0; i < n; i++) { + if (cells[i] == target) return true; + } + return false; +} + +SUITE(geodesicPolygonToCellsExperimental) { + sfGeoPolygon.geoloop = sfGeoLoop; + sfGeoPolygon.numHoles = 0; + + holeGeoPolygon.geoloop = sfGeoLoop; + holeGeoPolygon.numHoles = 1; + holeGeoPolygon.holes = &holeGeoLoop; + + nullHoleGeoPolygon.geoloop = sfGeoLoop; + nullHoleGeoPolygon.numHoles = 1; + nullHoleGeoPolygon.holes = &nullGeoLoop; + + pointHoleGeoPolygon.geoloop = sfGeoLoop; + pointHoleGeoPolygon.numHoles = 1; + pointHoleGeoPolygon.holes = &pointGeoLoop; + + lineHoleGeoPolygon.geoloop = sfGeoLoop; + lineHoleGeoPolygon.numHoles = 1; + lineHoleGeoPolygon.holes = &lineGeoLoop; + + emptyGeoPolygon.geoloop = emptyGeoLoop; + emptyGeoPolygon.numHoles = 0; + + invalidGeoPolygon.geoloop = invalidGeoLoop; + invalidGeoPolygon.numHoles = 0; + + invalid2GeoPolygon.geoloop = invalid2GeoLoop; + invalid2GeoPolygon.numHoles = 0; + + outOfBoundsVertGeoPolygon.geoloop = outOfBoundsVertGeoLoop; + outOfBoundsVertGeoPolygon.numHoles = 0; + + nullGeoPolygon.geoloop = nullGeoLoop; + nullGeoPolygon.numHoles = 0; + + pointGeoPolygon.geoloop = pointGeoLoop; + pointGeoPolygon.numHoles = 0; + + lineGeoPolygon.geoloop = lineGeoLoop; + lineGeoPolygon.numHoles = 0; + + wideWrapGeoPolygon.geoloop = wideWrapGeoLoop; + wideWrapGeoPolygon.numHoles = 0; + + tallSliceGeoPolygon.geoloop = tallSliceGeoLoop; + tallSliceGeoPolygon.numHoles = 0; + + reversedTriangleGeoPolygon.geoloop = reversedTriangleGeoLoop; + reversedTriangleGeoPolygon.numHoles = 0; + + transGeoPolygon.geoloop = transGeoLoop; + transGeoPolygon.numHoles = 0; + + TEST(edgePointRejected) { + static LatLng point[] = {{-0.0002458237579169511, 0.12401960784313724}}; + static GeoLoop pointLoop = {.numVerts = 1, .verts = point}; + static GeoPolygon pointPolygon; + pointPolygon.geoloop = pointLoop; + pointPolygon.numHoles = 0; + + int64_t size = 0; + H3Error err = geodesicFill(&pointPolygon, 1, CONTAINMENT_OVERLAPPING, + &size, NULL); + t_assert(err == E_DOMAIN, "single-vertex geodesic loops are rejected"); + } + + TEST(geodesicZeroSize) { + uint32_t flags = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(flags); + t_assert(H3_EXPORT(polygonToCellsExperimental)( + &sfGeoPolygon, 0, flags, 0, NULL) == E_MEMORY_BOUNDS, + "requires capacity for geodesic polyfill"); + } + + TEST(geodesicFullContainment) { + int64_t count = geodesicFillCount(&sfGeoPolygon, 5, CONTAINMENT_FULL); + t_assert(count == 2198, "matches expected geodesic full count"); + } + + TEST(geodesicOverlapping) { + int64_t count = + geodesicFillCount(&sfGeoPolygon, 5, CONTAINMENT_OVERLAPPING); + t_assert(count == 2469, "matches expected geodesic overlapping count"); + } + + TEST(geodesicCenter) { + int64_t fullCount = + geodesicFillCount(&sfGeoPolygon, 5, CONTAINMENT_FULL); + int64_t centerCount = + geodesicFillCount(&sfGeoPolygon, 5, CONTAINMENT_CENTER); + int64_t overlapCount = + geodesicFillCount(&sfGeoPolygon, 5, CONTAINMENT_OVERLAPPING); + t_assert(centerCount >= fullCount, + "center count should be >= geodesic full count"); + t_assert(centerCount <= overlapCount, + "center count should be <= geodesic overlapping count"); + } + + TEST(geodesicHoleFullContainment) { + int64_t count = geodesicFillCount(&holeGeoPolygon, 5, CONTAINMENT_FULL); + t_assert(count == 2123, "hole respected in full containment"); + } + + TEST(geodesicHoleOverlapping) { + int64_t count = + geodesicFillCount(&holeGeoPolygon, 5, CONTAINMENT_OVERLAPPING); + t_assert(count == 2436, "hole respected in overlapping containment"); + } + + TEST(geodesicHoleCenter) { + int64_t fullCount = + geodesicFillCount(&holeGeoPolygon, 5, CONTAINMENT_FULL); + int64_t centerCount = + geodesicFillCount(&holeGeoPolygon, 5, CONTAINMENT_CENTER); + int64_t overlapCount = + geodesicFillCount(&holeGeoPolygon, 5, CONTAINMENT_OVERLAPPING); + t_assert(centerCount >= fullCount, + "center count should be >= geodesic full count with holes"); + t_assert(centerCount <= overlapCount, + "center count should be <= geodesic overlap count with holes"); + } + + TEST(geodesicTransmeridianFull) { + int64_t count = + geodesicFillCount(&transGeoPolygon, 4, CONTAINMENT_FULL); + t_assert(count == 2246, "large transmeridian full count"); + } + + TEST(geodesicTransmeridianOverlapping) { + int64_t count = + geodesicFillCount(&transGeoPolygon, 4, CONTAINMENT_OVERLAPPING); + t_assert(count == 2554, "large transmeridian overlapping count"); + } + + TEST(geodesicTransmeridianCenter) { + int64_t fullCount = + geodesicFillCount(&transGeoPolygon, 4, CONTAINMENT_FULL); + int64_t centerCount = + geodesicFillCount(&transGeoPolygon, 4, CONTAINMENT_CENTER); + int64_t overlapCount = + geodesicFillCount(&transGeoPolygon, 4, CONTAINMENT_OVERLAPPING); + t_assert(centerCount >= fullCount, + "center count should be >= geodesic full count transmeridian"); + t_assert(centerCount <= overlapCount, + "center count should be <= geodesic overlapping count " + "transmeridian"); + } + + TEST(geodesicEmptyPolygon) { + int64_t count = + geodesicFillCount(&emptyGeoPolygon, 5, CONTAINMENT_FULL); + t_assert(count == 0, "empty polygon yields no cells"); + } + + TEST(geodesicNullPolygon) { + int64_t count = + geodesicFillCount(&nullGeoPolygon, 5, CONTAINMENT_OVERLAPPING); + t_assert(count == 0, "null polygon yields no cells"); + } + + TEST(geodesicPointPolygon) { + int64_t size = 0; + H3Error err = geodesicFill(&pointGeoPolygon, 5, CONTAINMENT_OVERLAPPING, + &size, NULL); + t_assert(err == E_DOMAIN, "single-vertex polygon rejected"); + } + + TEST(geodesicLinePolygon) { + int64_t size = 0; + H3Error err = geodesicFill(&lineGeoPolygon, 5, CONTAINMENT_OVERLAPPING, + &size, NULL); + t_assert(err == E_DOMAIN, "two-vertex polygon rejected"); + } + + TEST(geodesicInvalidLoopValues) { + int64_t size = 0; + H3Error err = + geodesicFill(&invalidGeoPolygon, 5, CONTAINMENT_FULL, &size, NULL); + t_assert(err == E_DOMAIN, "infinite verts rejected with E_DOMAIN"); + + err = + geodesicFill(&invalid2GeoPolygon, 5, CONTAINMENT_FULL, &size, NULL); + t_assert(err == E_DOMAIN, "NaN verts rejected with E_DOMAIN"); + } + + TEST(geodesicMismatchedHoles) { + // Test case: numHoles > 0 but holes pointer is NULL + GeoPolygon mismatchedPolygon = { + .geoloop = sfGeoLoop, + .numHoles = 1, // Says there's 1 hole + .holes = NULL // But holes pointer is NULL + }; + + int64_t size = 0; + H3Error err = + geodesicFill(&mismatchedPolygon, 5, CONTAINMENT_FULL, &size, NULL); + t_assert(err == E_DOMAIN, + "numHoles > 0 with NULL holes returns E_DOMAIN"); + } + + TEST(geodesicNullHolePolygon) { + // Test polygon with a hole that has 0 vertices + int64_t size = 0; + H3Error err = + geodesicFill(&nullHoleGeoPolygon, 5, CONTAINMENT_FULL, &size, NULL); + t_assert(err == E_DOMAIN, "hole with 0 vertices returns E_DOMAIN"); + } + + TEST(geodesicPointHolePolygon) { + int64_t size = 0; + H3Error err = geodesicFill(&pointHoleGeoPolygon, 5, + CONTAINMENT_OVERLAPPING, &size, NULL); + t_assert(err == E_DOMAIN, "single-vertex hole rejected"); + } + + TEST(geodesicLineHolePolygon) { + int64_t size = 0; + H3Error err = geodesicFill(&lineHoleGeoPolygon, 5, + CONTAINMENT_OVERLAPPING, &size, NULL); + t_assert(err == E_DOMAIN, "two-vertex hole rejected"); + } + + TEST(geodesicInvalidContainmentModes) { + int64_t size = 0; + H3Error err = geodesicFill(&sfGeoPolygon, 5, + CONTAINMENT_OVERLAPPING_BBOX, &size, NULL); + t_assert(err == E_OPTION_INVALID, + "overlapping bbox mode rejected for geodesic"); + err = geodesicFill(&sfGeoPolygon, 5, CONTAINMENT_INVALID, &size, NULL); + t_assert(err == E_OPTION_INVALID, + "invalid containment mode rejected for geodesic"); + err = geodesicFill(&sfGeoPolygon, 5, CONTAINMENT_INVALID + 1, &size, + NULL); + t_assert(err == E_OPTION_INVALID, + "out-of-range containment mode rejected for geodesic"); + } + + TEST(geodesicMaxSizeNoUnderAllocation) { + LatLng reproVerts[] = { + {H3_EXPORT(degsToRads)(20.0), H3_EXPORT(degsToRads)(-70.0)}, + {H3_EXPORT(degsToRads)(20.0), H3_EXPORT(degsToRads)(70.0)}, + {H3_EXPORT(degsToRads)(-5.0), H3_EXPORT(degsToRads)(0.0)}}; + GeoLoop reproGeoLoop = {.numVerts = 3, .verts = reproVerts}; + GeoPolygon reproGeoPolygon = {.geoloop = reproGeoLoop, .numHoles = 0}; + + uint32_t flags = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(flags); + + int64_t size = 0; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &reproGeoPolygon, 1, flags, &size)); + t_assert(size > 0, "max size estimate is non-zero"); + + H3Index *out = calloc(size, sizeof(H3Index)); + t_assert(out != NULL, "allocated output buffer"); + + H3Error err = H3_EXPORT(polygonToCellsExperimental)(&reproGeoPolygon, 1, + flags, size, out); + t_assert(err == E_SUCCESS, + "maxPolygonToCellsSizeExperimental must not under-allocate"); + + free(out); + } + + TEST(geodesicRejectsPolygonsLargerThanHemisphere) { + int64_t size = 0; + H3Error err = geodesicFill(&wideWrapGeoPolygon, 2, + CONTAINMENT_OVERLAPPING, &size, NULL); + t_assert(err == E_DOMAIN, + "very wide low-area loop that wraps globe is rejected"); + + err = geodesicFill(&tallSliceGeoPolygon, 2, CONTAINMENT_OVERLAPPING, + &size, NULL); + t_assert(err == E_DOMAIN, "very tall north-south loop is rejected"); + + err = geodesicFill(&reversedTriangleGeoPolygon, 2, + CONTAINMENT_OVERLAPPING, &size, NULL); + t_assert(err == E_DOMAIN, + "equatorial wrap loop that spans globe is rejected"); + } + + // Regression test: FULL and OVERLAPPING results must be identical for CCW + // and CW loop orientation. Previously, a too-tight epsilon in + // _geodesicEdgesCross caused coincident arcs (shared cell boundaries) to + // be misclassified, making OVERLAPPING sensitive to vertex order. Sweeps + // all res-0 (122 cells) and res-1 (842 cells) to catch any recurrence. + TEST(geodesicOrientationInvariantExhaustive) { + const int resolutions[] = {0, 1}; + const uint32_t modes[] = {CONTAINMENT_FULL, CONTAINMENT_OVERLAPPING}; + + for (int ri = 0; ri < 2; ri++) { + int res = resolutions[ri]; + for (IterCellsResolution it = iterInitRes(res); it.h; + iterStepRes(&it)) { + H3Index cell = it.h; + + CellBoundary bnd = {0}; + t_assertSuccess(H3_EXPORT(cellToBoundary)(cell, &bnd)); + + LatLng ccwVerts[MAX_CELL_BNDRY_VERTS]; + LatLng cwVerts[MAX_CELL_BNDRY_VERTS]; + for (int i = 0; i < bnd.numVerts; i++) { + ccwVerts[i] = bnd.verts[i]; + cwVerts[i] = bnd.verts[bnd.numVerts - 1 - i]; + } + + GeoLoop loopCCW = {.numVerts = bnd.numVerts, .verts = ccwVerts}; + GeoLoop loopCW = {.numVerts = bnd.numVerts, .verts = cwVerts}; + GeoPolygon polyCCW = {.geoloop = loopCCW, .numHoles = 0}; + GeoPolygon polyCW = {.geoloop = loopCW, .numHoles = 0}; + + for (int mi = 0; mi < 2; mi++) { + H3Index *outCCW = NULL, *outCW = NULL; + int64_t nCCW = + _geodesicFillSorted(&polyCCW, res, modes[mi], &outCCW); + int64_t nCW = + _geodesicFillSorted(&polyCW, res, modes[mi], &outCW); + + t_assert(nCCW >= 0 && nCW >= 0, "fill succeeded"); + t_assert(nCCW == nCW, + "CCW and CW must yield same count for each mode"); + + bool equal = true; + for (int64_t i = 0; i < nCCW && equal; i++) { + if (outCCW[i] != outCW[i]) equal = false; + } + t_assert( + equal, + "CCW and CW must yield identical sorted cell sets"); + + free(outCCW); + free(outCW); + } + } + } + } + + TEST(geodesicRejectsGreatCircleBoundaryLoop) { + LatLng equatorLineVerts[] = {{0.0, -60.0 * M_PI / 180.0}, + {0.0, -20.0 * M_PI / 180.0}, + {0.0, 20.0 * M_PI / 180.0}, + {0.0, 60.0 * M_PI / 180.0}}; + GeoLoop equatorLineLoop = {.numVerts = ARRAY_SIZE(equatorLineVerts), + .verts = equatorLineVerts}; + GeoPolygon equatorLinePolygon = { + .geoloop = equatorLineLoop, .numHoles = 0, .holes = NULL}; + + int64_t size = 0; + H3Error err = geodesicFill(&equatorLinePolygon, 2, + CONTAINMENT_OVERLAPPING, &size, NULL); + t_assert(err == E_DOMAIN, + "all-vertices-on-great-circle loop is rejected"); + } + + // Regression: A polygon far smaller than a cell contains that cell's center + TEST(geodesicFullRejectsPolygonSmallerThanCell) { + const int res = 5; + LatLng seed = {H3_EXPORT(degsToRads)(37.5), + H3_EXPORT(degsToRads)(-122.0)}; + H3Index cell; + t_assertSuccess(H3_EXPORT(latLngToCell)(&seed, res, &cell)); + + LatLng center; + t_assertSuccess(H3_EXPORT(cellToLatLng)(cell, ¢er)); + + // Tiny triangle centered on the cell center, well inside the cell + const double d = 0.0002; // radians, ~1.3 km + LatLng triVerts[] = {{center.lat + d, center.lng}, + {center.lat - d, center.lng + d}, + {center.lat - d, center.lng - d}}; + GeoLoop triLoop = {.numVerts = 3, .verts = triVerts}; + GeoPolygon triPolygon = {.geoloop = triLoop, .numHoles = 0}; + + int64_t overlapCount = + geodesicFillCount(&triPolygon, res, CONTAINMENT_OVERLAPPING); + t_assert(overlapCount >= 1, "tiny polygon overlaps its enclosing cell"); + + int64_t fullCount = + geodesicFillCount(&triPolygon, res, CONTAINMENT_FULL); + t_assert(fullCount == 0, + "a polygon smaller than a cell fully contains no cell"); + } + + // Regression: A small hole that does not contain the cell center + TEST(geodesicFullRespectsHoleInsideCell) { + const int res = 5; + LatLng seed = {H3_EXPORT(degsToRads)(37.5), + H3_EXPORT(degsToRads)(-117.5)}; + H3Index cell; + t_assertSuccess(H3_EXPORT(latLngToCell)(&seed, res, &cell)); + + LatLng center; + t_assertSuccess(H3_EXPORT(cellToLatLng)(cell, ¢er)); + + // A tiny hole offsetnorth of the cell center: entirely inside + // the cell, but not covering the center, and crossing no cell edge. + const double off = 0.0005; + const double hr = 0.0001; + LatLng holeVerts2[] = {{center.lat + off + hr, center.lng}, + {center.lat + off - hr, center.lng + hr}, + {center.lat + off - hr, center.lng - hr}}; + GeoLoop tinyHole = {.numVerts = 3, .verts = holeVerts2}; + GeoPolygon holed = { + .geoloop = sfGeoLoop, .numHoles = 1, .holes = &tinyHole}; + + H3Index *noHole = NULL; + H3Index *withHole = NULL; + int64_t nNo = + _geodesicFillSorted(&sfGeoPolygon, res, CONTAINMENT_FULL, &noHole); + int64_t nWith = + _geodesicFillSorted(&holed, res, CONTAINMENT_FULL, &withHole); + t_assert(nNo > 0 && nWith >= 0, "geodesic fills succeeded"); + + t_assert(_containsCell(noHole, nNo, cell), + "cell is fully contained when there is no hole"); + t_assert(!_containsCell(withHole, nWith, cell), + "cell is not fully contained once a hole sits inside it"); + + free(noHole); + free(withHole); + } +} diff --git a/src/apps/testapps/testGetIcosahedronFaces.c b/src/apps/testapps/testGetIcosahedronFaces.c index 73db7c0c68..8e42b0a95c 100644 --- a/src/apps/testapps/testGetIcosahedronFaces.c +++ b/src/apps/testapps/testGetIcosahedronFaces.c @@ -21,6 +21,7 @@ #include #include "baseCells.h" +#include "faceijk.h" #include "h3Index.h" #include "h3api.h" #include "test.h" diff --git a/src/apps/testapps/testH3IndexInternal.c b/src/apps/testapps/testH3IndexInternal.c index c7ed6c5cbe..a498614049 100644 --- a/src/apps/testapps/testH3IndexInternal.c +++ b/src/apps/testapps/testH3IndexInternal.c @@ -19,6 +19,7 @@ * usage: `testH3IndexInternal` */ +#include #include #include #include @@ -27,6 +28,7 @@ #include "h3Index.h" #include "test.h" #include "utility.h" +#include "vec3d.h" SUITE(h3IndexInternal) { TEST(faceIjkToH3ExtremeCoordinates) { @@ -51,4 +53,48 @@ SUITE(h3IndexInternal) { FaceIJK fijk2K = {8, {2, 0, 20}}; t_assert(_faceIjkToH3(&fijk2K, 2) == 0, "k out of bounds at res 2"); } + + TEST(vec3ToCell_invalidInputs) { + Vec3d v = {1.0, 0.0, 0.0}; + H3Index out; + + // Test invalid resolution (negative) + t_assert(vec3ToCell(&v, -1, &out) == E_RES_DOMAIN, + "negative resolution rejected"); + + // Test invalid resolution (too high) + t_assert(vec3ToCell(&v, MAX_H3_RES + 1, &out) == E_RES_DOMAIN, + "resolution above MAX_H3_RES rejected"); + + // Test non-finite coordinates + Vec3d nanVec = {NAN, 0.0, 0.0}; + t_assert(vec3ToCell(&nanVec, 5, &out) == E_DOMAIN, + "NaN x coordinate rejected"); + + Vec3d infVec = {INFINITY, 0.0, 0.0}; + t_assert(vec3ToCell(&infVec, 5, &out) == E_DOMAIN, + "Infinity x coordinate rejected"); + + Vec3d nanY = {0.0, NAN, 0.0}; + t_assert(vec3ToCell(&nanY, 5, &out) == E_DOMAIN, + "NaN y coordinate rejected"); + + Vec3d nanZ = {0.0, 0.0, NAN}; + t_assert(vec3ToCell(&nanZ, 5, &out) == E_DOMAIN, + "NaN z coordinate rejected"); + } + + TEST(cellToVec3_invalidCell) { + Vec3d v; + H3Index invalid = 0xFFFFFFFFFFFFFFFF; + t_assert(cellToVec3(invalid, &v) != E_SUCCESS, + "invalid H3 index rejected in cellToVec3"); + } + + TEST(cellToGeodesicBoundary_invalidCell) { + GeodesicCellBoundary cb; + H3Index invalid = 0xFFFFFFFFFFFFFFFF; + t_assert(cellToGeodesicBoundary(invalid, &cb) != E_SUCCESS, + "invalid H3 index rejected in cellToGeodesicBoundary"); + } } diff --git a/src/apps/testapps/testH3Memory.c b/src/apps/testapps/testH3Memory.c index ce934614fe..8e6239c850 100644 --- a/src/apps/testapps/testH3Memory.c +++ b/src/apps/testapps/testH3Memory.c @@ -13,8 +13,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -/** @file - * @brief tests main H3 core library memory management. +/** @file testH3Memory.c + * @brief Exercises memory allocation hooks in the H3 core library. * * usage: `testH3Memory` */ @@ -23,6 +23,8 @@ #include #include "cellsToMultiPoly.h" +#include "geodesicIterator.h" +#include "geodesicPolygonInternal.h" #include "h3Index.h" #include "h3api.h" #include "latLng.h" @@ -551,4 +553,259 @@ SUITE(h3Memory) { "Should fail with memory error before success"); } } + + TEST(polygonToCellsExperimentalGeodesic) { + sfGeoPolygon.geoloop = sfGeoLoop; + sfGeoPolygon.numHoles = 0; + + uint32_t flags = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(flags); + + // Test maxPolygonToCellsSizeExperimental with allocation failure + resetMemoryCounters(0); + failAlloc = true; + int64_t numHexagons; + H3Error err = H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 5, flags, &numHexagons); + t_assert(err == E_MEMORY_ALLOC, + "geodesic maxPolygonToCellsSizeExperimental failed on alloc"); + t_assert(actualAllocCalls == 1, + "maxPolygonToCellsSizeExperimental alloc called once"); + t_assert(actualFreeCalls == 0, + "maxPolygonToCellsSizeExperimental free not called"); + + // Test maxPolygonToCellsSizeExperimental with allocation success. + // Geodesic sizing may perform multiple allocations internally. + resetMemoryCounters(0); + failAlloc = false; + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 5, flags, &numHexagons)); + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + + // Test polygonToCellsExperimental with allocation failure + resetMemoryCounters(0); + failAlloc = true; + err = H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 5, flags, + numHexagons, hexagons); + t_assert(err == E_MEMORY_ALLOC, + "geodesic polygonToCellsExperimental failed (1)"); + t_assert(actualAllocCalls == 1, "alloc called once"); + t_assert(actualFreeCalls == 0, "free not called"); + + resetMemoryCounters(3); + err = H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 5, flags, + numHexagons, hexagons); + t_assert(err == E_SUCCESS, + "geodesic polygonToCellsExperimental succeeded (1)"); + t_assert(actualAllocCalls == 3, "alloc called three times"); + t_assert(actualFreeCalls == 3, "free called three times"); + + int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons); + t_assert(actualNumIndexes == 3, + "got expected geodesic polygonToCellsExperimental size"); + free(hexagons); + } + + TEST(geodesicPolygonCreateAllocationFailures) { + LatLng outerVerts[] = {{0.659966917655, -2.1364398519396}, + {0.6595011102219, -2.1359434279405}, + {0.6583348114025, -2.1354884206045}}; + GeoLoop outer = {.numVerts = 3, .verts = outerVerts}; + + LatLng hole0Verts[] = { + {0.6598, -2.1362}, {0.6597, -2.1360}, {0.6596, -2.1362}}; + LatLng hole1Verts[] = { + {0.65955, -2.13630}, {0.65950, -2.13615}, {0.65945, -2.13630}}; + GeoLoop validHoles[] = {{.numVerts = 3, .verts = hole0Verts}, + {.numVerts = 3, .verts = hole1Verts}}; + GeoLoop secondHoleInvalid[] = {{.numVerts = 3, .verts = hole0Verts}, + {.numVerts = 0, .verts = NULL}}; + + GeoPolygon noHoles = {.geoloop = outer, .numHoles = 0, .holes = NULL}; + GeoPolygon oneHole = { + .geoloop = outer, .numHoles = 1, .holes = validHoles}; + GeoPolygon twoHolesInvalid = { + .geoloop = outer, .numHoles = 2, .holes = secondHoleInvalid}; + GeoPolygon twoHolesValid = { + .geoloop = outer, .numHoles = 2, .holes = validHoles}; + + GeodesicPolygon *poly = NULL; + + // Fail allocating the GeodesicPolygon container. + resetMemoryCounters(0); + failAlloc = true; + t_assert(geodesicPolygonCreate(&noHoles, &poly) == E_MEMORY_ALLOC, + "geodesicPolygonCreate fails with first allocation failure"); + t_assert(actualAllocCalls == 1, "first allocation attempted once"); + t_assert(actualFreeCalls == 0, "no free needed when first alloc fails"); + + // Fail allocating the outer loop edges. + resetMemoryCounters(1); + t_assert(geodesicPolygonCreate(&noHoles, &poly) == E_MEMORY_ALLOC, + "geodesicPolygonCreate fails creating outer loop"); + t_assert(actualAllocCalls == 2, + "outer loop failure reached the second allocation"); + t_assert(actualFreeCalls == 1, "container freed on outer loop failure"); + + // Fail allocating the holes array. + resetMemoryCounters(2); + t_assert(geodesicPolygonCreate(&oneHole, &poly) == E_MEMORY_ALLOC, + "geodesicPolygonCreate fails allocating holes array"); + t_assert(actualAllocCalls == 3, + "holes-array failure reached the third allocation"); + t_assert(actualFreeCalls == 2, "outer loop and container freed"); + + // Invalid second hole cleans up first hole and container state. + resetMemoryCounters(0); + t_assert(geodesicPolygonCreate(&twoHolesInvalid, &poly) == E_DOMAIN, + "invalid second hole triggers cleanup path"); + t_assert(actualAllocCalls == 4, "allocations reached first hole edges"); + t_assert(actualFreeCalls == 4, "all intermediate allocations freed"); + + // Fail allocating the second hole edges after first hole succeeds. + resetMemoryCounters(4); + t_assert(geodesicPolygonCreate(&twoHolesValid, &poly) == E_MEMORY_ALLOC, + "second hole allocation failure triggers deep cleanup"); + t_assert(actualAllocCalls == 5, + "second-hole failure reached fifth allocation"); + t_assert(actualFreeCalls == 4, + "first hole, holes array, outer loop, and container freed"); + } + + TEST(geodesicIteratorStepDefensivePaths) { + uint32_t geodesicFlags = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(geodesicFlags); + + LatLng pointVerts[] = {{0.659966917655, -2.1364398519396}}; + GeoPolygon pointLoop = {.geoloop = {.numVerts = 1, .verts = pointVerts}, + .numHoles = 0, + .holes = NULL}; + IterCellsPolygonCompact pointLoopIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &pointLoop, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + geodesicIteratorStep(&pointLoopIter, pointLoopIter.cell); + t_assert(pointLoopIter.error == E_DOMAIN, + "single-vertex geodesic iterator polygon rejected"); + t_assert(pointLoopIter.cell == H3_NULL, "iterator moved to exhausted"); + + LatLng lineVerts[] = {{0.659966917655, -2.1364398519396}, + {0.6595011102219, -2.1359434279405}}; + GeoPolygon lineLoop = {.geoloop = {.numVerts = 2, .verts = lineVerts}, + .numHoles = 0, + .holes = NULL}; + IterCellsPolygonCompact lineLoopIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &lineLoop, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + geodesicIteratorStep(&lineLoopIter, lineLoopIter.cell); + t_assert(lineLoopIter.error == E_DOMAIN, + "two-vertex geodesic iterator polygon rejected"); + t_assert(lineLoopIter.cell == H3_NULL, "iterator moved to exhausted"); + + GeoPolygon zeroLoop = {.geoloop = {.numVerts = 0, .verts = NULL}, + .numHoles = 0, + .holes = NULL}; + IterCellsPolygonCompact zeroLoopIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &zeroLoop, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + geodesicIteratorStep(&zeroLoopIter, zeroLoopIter.cell); + t_assert(zeroLoopIter.error == E_DOMAIN, + "zero-vertex geodesic iterator polygon rejected"); + t_assert(zeroLoopIter.cell == H3_NULL, "iterator moved to exhausted"); + + LatLng outerVerts[] = {{0.659966917655, -2.1364398519396}, + {0.6595011102219, -2.1359434279405}, + {0.6583348114025, -2.1354884206045}}; + GeoPolygon nullHoles = {.geoloop = {.numVerts = 3, .verts = outerVerts}, + .numHoles = 1, + .holes = NULL}; + IterCellsPolygonCompact nullHolesIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &nullHoles, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + geodesicIteratorStep(&nullHolesIter, nullHolesIter.cell); + t_assert(nullHolesIter.error == E_DOMAIN, + "holes count with NULL holes pointer rejected"); + t_assert(nullHolesIter.cell == H3_NULL, "iterator moved to exhausted"); + + GeoLoop pointHole = {.numVerts = 1, .verts = pointVerts}; + GeoPolygon pointHolePoly = { + .geoloop = {.numVerts = 3, .verts = outerVerts}, + .numHoles = 1, + .holes = &pointHole}; + IterCellsPolygonCompact pointHoleIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &pointHolePoly, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + geodesicIteratorStep(&pointHoleIter, pointHoleIter.cell); + t_assert(pointHoleIter.error == E_DOMAIN, + "single-vertex geodesic iterator hole rejected"); + t_assert(pointHoleIter.cell == H3_NULL, "iterator moved to exhausted"); + + GeoLoop lineHole = {.numVerts = 2, .verts = lineVerts}; + GeoPolygon lineHolePoly = { + .geoloop = {.numVerts = 3, .verts = outerVerts}, + .numHoles = 1, + .holes = &lineHole}; + IterCellsPolygonCompact lineHoleIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &lineHolePoly, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + geodesicIteratorStep(&lineHoleIter, lineHoleIter.cell); + t_assert(lineHoleIter.error == E_DOMAIN, + "two-vertex geodesic iterator hole rejected"); + t_assert(lineHoleIter.cell == H3_NULL, "iterator moved to exhausted"); + + GeoPolygon valid = {.geoloop = {.numVerts = 3, .verts = outerVerts}, + .numHoles = 0, + .holes = NULL}; + IterCellsPolygonCompact allocFailIter = {.cell = 0x8001fffffffffff, + .error = E_SUCCESS, + ._res = 0, + ._flags = geodesicFlags, + ._polygon = &valid, + ._bboxes = NULL, + ._started = true, + .geodesicPoly = NULL}; + resetMemoryCounters(0); + failAlloc = true; + geodesicIteratorStep(&allocFailIter, allocFailIter.cell); + t_assert(allocFailIter.error == E_MEMORY_ALLOC, + "geodesic iterator reports allocation failure"); + t_assert(allocFailIter.cell == H3_NULL, "iterator moved to exhausted"); + t_assert(actualAllocCalls == 1, + "attempted geodesic polygon allocation"); + t_assert(actualFreeCalls == 0, + "no state to free when geodesic allocation fails"); + + geodesicIteratorDestroyState(NULL); + geodesicIteratorDestroyState(&allocFailIter); + t_assert(true, "destroy state accepts null and empty iterator"); + } } diff --git a/src/apps/testapps/testPolygonToCells.c b/src/apps/testapps/testPolygonToCells.c index f8724bac91..0fb0eed36d 100644 --- a/src/apps/testapps/testPolygonToCells.c +++ b/src/apps/testapps/testPolygonToCells.c @@ -14,6 +14,10 @@ * limitations under the License. */ +/** @file testPolygonToCells.c + * @brief Tests the legacy polygon-to-cells traversal algorithm. + */ + #include #include #include diff --git a/src/apps/testapps/testPolygonToCellsExperimental.c b/src/apps/testapps/testPolygonToCellsExperimental.c index 371629d4e6..c4a57359a0 100644 --- a/src/apps/testapps/testPolygonToCellsExperimental.c +++ b/src/apps/testapps/testPolygonToCellsExperimental.c @@ -14,6 +14,10 @@ * limitations under the License. */ +/** @file testPolygonToCellsExperimental.c + * @brief Tests the experimental polygon-to-cells traversal algorithm. + */ + #include #include #include @@ -1090,23 +1094,99 @@ SUITE(polygonToCells) { TEST(invalidFlags) { int64_t numHexagons; - for (uint32_t flags = CONTAINMENT_INVALID; flags <= 32; flags++) { + + // Test invalid containment modes (>= CONTAINMENT_INVALID) for + // non-geodesic + for (uint32_t mode = CONTAINMENT_INVALID; + mode <= CONTAINMENT_INVALID + 5; mode++) { t_assert( H3_EXPORT(maxPolygonToCellsSizeExperimental)( - &sfGeoPolygon, 9, flags, &numHexagons) == E_OPTION_INVALID, - "Flags other than polyfill modes are invalid for " - "maxPolygonToCellsSizeExperimental"); + &sfGeoPolygon, 9, mode, &numHexagons) == E_OPTION_INVALID, + "Invalid containment modes should be rejected"); } + + // Test invalid containment modes for geodesic + uint32_t invalidGeodesicFlags[] = { + FLAG_GEODESIC_MASK | CONTAINMENT_OVERLAPPING_BBOX, + FLAG_GEODESIC_MASK | CONTAINMENT_INVALID}; + for (int i = 0; i < 2; i++) { + t_assert(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, invalidGeodesicFlags[i], + &numHexagons) == E_OPTION_INVALID, + "Invalid geodesic containment modes should be rejected"); + } + + // Test invalid flag bits (bits outside valid masks) + t_assert(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, 0x100, &numHexagons) == E_OPTION_INVALID, + "Invalid flag bits should be rejected"); + + // Test valid flags work t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons)); + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, FLAG_GEODESIC_MASK | CONTAINMENT_FULL, + &numHexagons)); + + // Test polygonToCellsExperimental with same invalid flags H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); - for (uint32_t flags = CONTAINMENT_INVALID; flags <= 32; flags++) { - t_assert(H3_EXPORT(polygonToCellsExperimental)( - &sfGeoPolygon, 9, flags, numHexagons, hexagons) == - E_OPTION_INVALID, - "Flags other than polyfill modes are invalid for " - "polygonToCellsExperimental"); + t_assert(H3_EXPORT(polygonToCellsExperimental)( + &sfGeoPolygon, 9, CONTAINMENT_INVALID, numHexagons, + hexagons) == E_OPTION_INVALID, + "Invalid flags should be rejected"); + t_assert(H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0x100, + numHexagons, hexagons) == + E_OPTION_INVALID, + "Invalid flag bits should be rejected"); + free(hexagons); + } + + TEST(invalidFlagsExperimental) { + int64_t numHexagons; + + // Test invalid containment modes (>= CONTAINMENT_INVALID) for + // non-geodesic + for (uint32_t mode = CONTAINMENT_INVALID; + mode <= CONTAINMENT_INVALID + 5; mode++) { + t_assert( + H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, mode, &numHexagons) == E_OPTION_INVALID, + "Invalid containment modes should be rejected"); } + + // Test invalid containment modes for geodesic + uint32_t invalidGeodesicFlags[] = { + FLAG_GEODESIC_MASK | CONTAINMENT_OVERLAPPING_BBOX, + FLAG_GEODESIC_MASK | CONTAINMENT_INVALID}; + for (int i = 0; i < 2; i++) { + t_assert(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, invalidGeodesicFlags[i], + &numHexagons) == E_OPTION_INVALID, + "Invalid geodesic containment modes should be rejected"); + } + + // Test invalid flag bits (bits outside valid masks) + t_assert(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, 0x100, &numHexagons) == E_OPTION_INVALID, + "Invalid flag bits should be rejected"); + + // Test valid flags work + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons)); + t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)( + &sfGeoPolygon, 9, FLAG_GEODESIC_MASK | CONTAINMENT_FULL, + &numHexagons)); + + // Test polygonToCellsExperimental with same invalid flags + H3Index *hexagons = calloc(numHexagons, sizeof(H3Index)); + t_assert(H3_EXPORT(polygonToCellsExperimental)( + &sfGeoPolygon, 9, CONTAINMENT_INVALID, numHexagons, + hexagons) == E_OPTION_INVALID, + "Invalid flags should be rejected"); + t_assert(H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0x100, + numHexagons, hexagons) == + E_OPTION_INVALID, + "Invalid flag bits should be rejected"); free(hexagons); } diff --git a/src/apps/testapps/testSphereCap.c b/src/apps/testapps/testSphereCap.c new file mode 100644 index 0000000000..20e27e8321 --- /dev/null +++ b/src/apps/testapps/testSphereCap.c @@ -0,0 +1,199 @@ +/* + * Copyright 2016-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. + */ + +/** @file testSphereCap.c + * @brief Tests the internal SphereCap helpers used by geodesic algorithms + */ + +#include +#include + +#include "bbox.h" +#include "constants.h" +#include "h3Index.h" +#include "latLng.h" +#include "polyfill.h" +#include "sphereCapTables.h" +#include "test.h" +#include "vec3d.h" + +// If update of the precomputed values if needed for any reason this +// function can be used to re-record them +static void _printPrecomputedCosRadius(void) { + printf("Printing \n"); + for (int res = 0; res <= MAX_H3_RES; res++) { + printf("res=%2d: %a\n", res, + cos(MAX_EDGE_LENGTH_RADS[res] * H3_SPHERE_CAP_SCALE_FACTOR)); + } +} + +/** + * Verify that PRECOMPUTED_COS_RADIUS matches cos(MAX_EDGE_LENGTH_RADS[res] * + * H3_SPHERE_CAP_SCALE_FACTOR) for every resolution. + */ +static bool _verifyCosRadiusValues(void) { + bool pass = true; + + for (int res = 0; res <= MAX_H3_RES; res++) { + const double radius = + MAX_EDGE_LENGTH_RADS[res] * H3_SPHERE_CAP_SCALE_FACTOR; + const double expected = cos(radius); + const double actual = PRECOMPUTED_COS_RADIUS[res]; + const double diff = fabs(expected - actual); + + if (diff > EPSILON) { + printf( + "Cos radius mismatch at res=%d: expected=%a actual=%a " + "diff=%0.3e\n", + res, expected, actual, diff); + pass = false; + } + } + + if (!pass) { + _printPrecomputedCosRadius(); + } + + return pass; +} + +static bool _cellWithinCap(const H3Index cell, const SphereCap *cap) { + CellBoundary boundary; + if (H3_EXPORT(cellToBoundary)(cell, &boundary) != E_SUCCESS) { + return false; + } + + for (int i = 0; i < boundary.numVerts; i++) { + Vec3d vertex; + vertex = latLngToVec3(boundary.verts[i]); + double dot = vec3Dot(vertex, cap->center); + if (dot + EPSILON < cap->cosRadius) { + printf( + "Cell %llx vertex %d outside cap (dot=%0.6f, thresh=%0.6f)\n", + (unsigned long long)cell, i, dot, cap->cosRadius); + return false; + } + + Vec3d next; + next = latLngToVec3(boundary.verts[(i + 1) % boundary.numVerts]); + Vec3d midpoint = {vertex.x + next.x, vertex.y + next.y, + vertex.z + next.z}; + double mag = vec3Norm(midpoint); + if (mag > EPSILON) { + midpoint.x /= mag; + midpoint.y /= mag; + midpoint.z /= mag; + double midDot = vec3Dot(midpoint, cap->center); + if (midDot + EPSILON < cap->cosRadius) { + printf( + "Cell %llx edge %d midpoint outside cap (dot=%0.6f, " + "thresh=%0.6f)\n", + (unsigned long long)cell, i, midDot, cap->cosRadius); + return false; + } + } + } + + return true; +} + +static bool _testCellRecursive(const H3Index cell, const int minTestRes, + const int maxTestRes) { + SphereCap cap; + if (cellToSphereCap(cell, &cap) != E_SUCCESS) { + return false; + } + + if (!_cellWithinCap(cell, &cap)) { + return false; + } + + int cellRes = H3_GET_RESOLUTION(cell); + if (cellRes >= maxTestRes) { + return true; + } + + int64_t childCount = 0; + if (H3_EXPORT(cellToChildrenSize)(cell, cellRes + 1, &childCount) != + E_SUCCESS || + childCount == 0) { + return false; + } + + H3Index *children = calloc(childCount, sizeof(H3Index)); + if (!children) { + return false; + } + + H3Error err = H3_EXPORT(cellToChildren)(cell, cellRes + 1, children); + if (err != E_SUCCESS) { + free(children); + return false; + } + + bool pass = true; + for (int64_t i = 0; i < childCount; i++) { + if (cellRes >= minTestRes) { + if (!_cellWithinCap(children[i], &cap)) { + pass = false; + break; + } + } + if (!_testCellRecursive(children[i], minTestRes, maxTestRes)) { + pass = false; + break; + } + } + + free(children); + return pass; +} + +SUITE(SphereCap) { + TEST(cellToSphereCapOutputs) { + const LatLng sf = {0.659966917655, -2.1364398519394}; + H3Index cell; + t_assertSuccess(H3_EXPORT(latLngToCell)(&sf, 3, &cell)); + + SphereCap cap; + t_assertSuccess(cellToSphereCap(cell, &cap)); + t_assert(fabs(cap.cosRadius - 0x1.fff685c0bd115p-1) < 1e-15, + "Matches precomputed cosine bound at resolution 3"); + t_assert(fabs(vec3Norm(cap.center) - 1.0) < 1e-12, + "Center vector is normalized"); + } + + TEST(testCap) { + const int testRes = 0; + const int maxTestRes = 4; + + if (testRes < 0 || testRes > MAX_H3_RES || maxTestRes < testRes || + maxTestRes > MAX_H3_RES) { + t_assert(false, "testCaps invalid resolution parameters"); + } + + if (!_verifyCosRadiusValues()) { + t_assert(false, "Precomputed cos radii failed verification"); + } + + for (int baseCell = 0; baseCell < NUM_BASE_CELLS; baseCell++) { + H3Index baseIndex = baseCellNumToCell(baseCell); + if (!_testCellRecursive(baseIndex, testRes, maxTestRes)) { + t_assert(false, "Bounding cap test failed"); + } + } + } +} diff --git a/src/apps/testapps/testVec3d.c b/src/apps/testapps/testVec3d.c new file mode 100644 index 0000000000..e2aef5ffde --- /dev/null +++ b/src/apps/testapps/testVec3d.c @@ -0,0 +1,75 @@ +/* + * Copyright 2026 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 testVec3d.c + * @brief Tests the vec3d helpers used by the geodesic polyfill path. + */ + +#include + +#include "constants.h" +#include "test.h" +#include "vec3d.h" + +SUITE(Vec3d) { + TEST(dotProduct) { + Vec3d a = {.x = 1.0, .y = 0.0, .z = 0.0}; + Vec3d b = {.x = -1.0, .y = 0.0, .z = 0.0}; + t_assert(vec3Dot(a, b) == -1.0, "dot product matches expected value"); + } + + TEST(crossProductOrthogonality) { + Vec3d i = {.x = 1.0, .y = 0.0, .z = 0.0}; + Vec3d j = {.x = 0.0, .y = 1.0, .z = 0.0}; + Vec3d k; + k = vec3Cross(i, j); + t_assert(fabs(k.x - 0.0) < EPSILON, "x component zero"); + t_assert(fabs(k.y - 0.0) < EPSILON, "y component zero"); + t_assert(fabs(k.z - 1.0) < EPSILON, "z component one"); + t_assert(fabs(vec3Dot(k, i)) < EPSILON, "cross is orthogonal to i"); + t_assert(fabs(vec3Dot(k, j)) < EPSILON, "cross is orthogonal to j"); + } + + TEST(normalizeAndMagnitude) { + Vec3d v = {.x = 3.0, .y = -4.0, .z = 12.0}; + double magSq = vec3NormSq(v); + t_assert(fabs(magSq - 169.0) < EPSILON, "magnitude squared matches"); + t_assert(fabs(vec3Norm(v) - 13.0) < EPSILON, "magnitude matches"); + + vec3Normalize(&v); + t_assert(fabs(vec3Norm(v) - 1.0) < 1e-12, "normalized vector is unit"); + + Vec3d zero = {.x = 0.0, .y = 0.0, .z = 0.0}; + vec3Normalize(&zero); + t_assert(zero.x == 0.0 && zero.y == 0.0 && zero.z == 0.0, + "zero vector remains unchanged when normalizing"); + } + + TEST(distance) { + Vec3d a = {.x = 0.0, .y = 0.0, .z = 0.0}; + Vec3d b = {.x = 1.0, .y = 2.0, .z = 2.0}; + t_assert(fabs(vec3DistSq(a, b) - 9.0) < EPSILON, + "distance squared matches"); + } + + TEST(latLngConversionConsistency) { + LatLng geo = {.lat = 0.5, .lng = -1.3}; + Vec3d v = latLngToVec3(geo); + + t_assert(fabs(vec3Norm(v) - 1.0) < 1e-12, + "converted vector lives on the unit sphere"); + } +} diff --git a/src/h3lib/include/bbox.h b/src/h3lib/include/bbox.h index 1c56982f4c..3f242b4de5 100644 --- a/src/h3lib/include/bbox.h +++ b/src/h3lib/include/bbox.h @@ -24,6 +24,7 @@ #include "h3api.h" #include "latLng.h" +#include "vec3d.h" /** @struct BBox * @brief Geographic bounding box with coordinates defined in radians @@ -35,6 +36,23 @@ typedef struct { double west; ///< west longitude } BBox; +/** @struct AABB + * @brief Axis-aligned bounding box expressed in Cartesian space on the unit + * sphere. + */ +typedef struct { + Vec3d min; ///< Minimum corner of the box + Vec3d max; ///< Maximum corner of the box +} AABB; + +/** @struct SphereCap + * @brief Bounding cap described by its center vector and cosine radius. + */ +typedef struct { + Vec3d center; ///< Unit vector pointing to the cap center + double cosRadius; ///< Cosine of the angular radius of the cap +} SphereCap; + double bboxWidthRads(const BBox *bbox); double bboxHeightRads(const BBox *bbox); bool bboxIsTransmeridian(const BBox *bbox); @@ -52,4 +70,30 @@ void bboxNormalization(const BBox *a, const BBox *b, LongitudeNormalization *aNormalization, LongitudeNormalization *bNormalization); +/** + * Expand an AABB with extrema from the great-circle arc connecting two points. + * + * @param aabb Axis-aligned bounding box to expand. + * @param v1 First endpoint of the arc. + * @param v2 Second endpoint of the arc. + * @param n Normal vector of the great circle defined by the arc. + */ +void aabbUpdateWithArcExtrema(AABB *aabb, const Vec3d *v1, const Vec3d *v2, + const Vec3d *n); + +/** Reset an AABB to an empty inverted state (min > max on every axis). */ +void aabbEmptyInverted(AABB *box); + +/** Expand an AABB with a single Cartesian point. */ +void aabbUpdateWithVec3d(AABB *aabb, const Vec3d *v); + +/** + * Compute a bounding sphere cap for a cell index. + * + * @param cell Cell to bound. + * @param out Output bounding cap. + * @return `E_SUCCESS` on success, or another error code on failure. + */ +H3Error cellToSphereCap(H3Index cell, SphereCap *out); + #endif diff --git a/src/h3lib/include/faceijk.h b/src/h3lib/include/faceijk.h index 77ef642bc9..5e24645960 100644 --- a/src/h3lib/include/faceijk.h +++ b/src/h3lib/include/faceijk.h @@ -25,6 +25,7 @@ #define FACEIJK_H #include "coordijk.h" +#include "geodesicCellBoundary.h" #include "latLng.h" #include "vec2d.h" #include "vec3d.h" @@ -75,8 +76,12 @@ void _vec3ToFaceIjk(Vec3d p, int res, FaceIJK *h); void _faceIjkToVec3(const FaceIJK *h, int res, Vec3d *g); void _faceIjkToCellBoundary(const FaceIJK *h, int res, int start, int length, CellBoundary *g); +void _faceIjkToCellBoundaryGeodesic(const FaceIJK *h, int res, int start, + int length, GeodesicCellBoundary *g); void _faceIjkPentToCellBoundary(const FaceIJK *h, int res, int start, int length, CellBoundary *g); +void _faceIjkPentToCellBoundaryGeodesic(const FaceIJK *h, int res, int start, + int length, GeodesicCellBoundary *g); void _faceIjkToVerts(FaceIJK *fijk, int *res, FaceIJK *fijkVerts); void _faceIjkPentToVerts(FaceIJK *fijk, int *res, FaceIJK *fijkVerts); Overage _adjustOverageClassII(FaceIJK *fijk, int res, int pentLeading4, diff --git a/src/h3lib/include/geodesicCellBoundary.h b/src/h3lib/include/geodesicCellBoundary.h new file mode 100644 index 0000000000..6bf6458536 --- /dev/null +++ b/src/h3lib/include/geodesicCellBoundary.h @@ -0,0 +1,35 @@ +/* + * Copyright 2026 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 geodesicCellBoundary.h + * @brief Definition of geodesic cell boundary + */ + +#ifndef GEODESIC_CELL_BOUNDARY_H +#define GEODESIC_CELL_BOUNDARY_H + +#include "constants.h" +#include "h3api.h" // For MAX_CELL_BNDRY_VERTS +#include "vec3d.h" + +/** @struct GeodesicCellBoundary + @brief Geodesic cell boundary expressed in Cartesian space +*/ +typedef struct { + int numVerts; ///< number of vertices + Vec3d verts[MAX_CELL_BNDRY_VERTS]; ///< vertices in ccw order +} GeodesicCellBoundary; + +#endif diff --git a/src/h3lib/include/geodesicIterator.h b/src/h3lib/include/geodesicIterator.h new file mode 100644 index 0000000000..17c818ec2d --- /dev/null +++ b/src/h3lib/include/geodesicIterator.h @@ -0,0 +1,42 @@ +/* + * Copyright 2026 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 geodesicIterator.h + * @brief Declares helpers for the geodesic polygon iterator state machine. + */ + +#ifndef GEODESIC_ITERATOR_H +#define GEODESIC_ITERATOR_H + +#include "h3api.h" +#include "polyfillIterator.h" + +/** + * Advance the geodesic polyfill iterator state for the provided starting cell. + * + * @param iter Iterator that owns the geodesic acceleration structure. + * @param cell The cell to continue iterating from. + */ +void geodesicIteratorStep(IterCellsPolygonCompact *iter, H3Index cell); + +/** + * Release any geodesic polygon state allocated for the iterator. + * + * @param iter Iterator that may own geodesic acceleration state. + */ +void geodesicIteratorDestroyState(IterCellsPolygonCompact *iter); + +#endif diff --git a/src/h3lib/include/geodesicPolygonInternal.h b/src/h3lib/include/geodesicPolygonInternal.h new file mode 100644 index 0000000000..d29f2bfd76 --- /dev/null +++ b/src/h3lib/include/geodesicPolygonInternal.h @@ -0,0 +1,108 @@ +/* + * Copyright 2026 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 geodesicPolygonInternal.h + * @brief Internal data structures for geodesic polygon acceleration. + */ + +#ifndef GEODESIC_POLYGON_INTERNAL_H +#define GEODESIC_POLYGON_INTERNAL_H + +#include + +#include "bbox.h" +#include "geodesicCellBoundary.h" +#include "polygon.h" +#include "vec3d.h" + +/** @struct GeodesicEdge + * @brief Cached data describing a polygon edge on the unit sphere. + */ +typedef struct { + Vec3d vert; ///< Vertex position on the unit sphere + Vec3d edgeCross; ///< Great-circle normal for the edge starting at vert + double edgeDot; ///< Dot product between consecutive vertices + AABB aabb; ///< Bounding box covering the edge arc +} GeodesicEdge; + +/** @struct GeodesicLoop + * @brief Polygon loop represented in geodesic space. + */ +typedef struct { + int numEdges; ///< Number of edges in the loop + GeodesicEdge *edges; ///< Edge data for the loop + Vec3d centroid; ///< Approximate centroid used for hemisphere tests + Vec3d hemispherePole; ///< Pole for a containing hemisphere (if available) + bool hasHemispherePole; ///< Whether hemispherePole is initialized +} GeodesicLoop; + +/** @struct GeodesicPolygon + * @brief Geodesic acceleration structure for polygon containment tests. + */ +typedef struct { + GeodesicLoop geoloop; ///< Exterior boundary in geodesic space + AABB aabb; ///< Bounding box covering the entire polygon + int numHoles; ///< Number of interior loops + GeodesicLoop *holes; ///< Optional interior loops in geodesic space +} GeodesicPolygon; + +/** + * Build a geodesic acceleration structure from a geographic polygon. + * + * @param polygon Source polygon expressed in lat/lng. + * @param out Output pointer to the allocated geodesic polygon. + * @return E_SUCCESS on success, or an error code (e.g. E_DOMAIN for + * non-finite vertices, E_MEMORY_ALLOC for allocation failure). + */ +H3Error geodesicPolygonCreate(const GeoPolygon *polygon, GeodesicPolygon **out); + +/** Destroy a geodesic polygon created by `geodesicPolygonCreate`. */ +void geodesicPolygonDestroy(GeodesicPolygon *polygon); + +/** + * Check whether a cell bounding cap intersects the polygon's outer hull. + * + * @param polygon Geodesic polygon to test. + * @param cap Bounding cap derived from a candidate cell. + * @return `true` when the cap may intersect the polygon. + */ +bool geodesicPolygonCapIntersects(const GeodesicPolygon *polygon, + const SphereCap *cap); + +/** + * Test whether a cell boundary intersects any polygon loop. + * + * @param polygon Geodesic polygon to test. + * @param boundary Geodesic boundary of the candidate cell. + * @param cap Bounding cap for the candidate cell. + * @return `true` if the cell boundary crosses the polygon. + */ +bool geodesicPolygonBoundaryIntersects(const GeodesicPolygon *polygon, + const GeodesicCellBoundary *boundary, + const SphereCap *cap); + +/** + * Determine whether a point on the unit sphere lies inside the polygon. + * + * @param polygon Geodesic polygon to test. + * @param point Unit vector for the point of interest. + * @return `true` when the point is inside the polygon exterior and outside all + * holes. + */ +bool geodesicPolygonContainsPoint(const GeodesicPolygon *polygon, + const Vec3d *point); + +#endif diff --git a/src/h3lib/include/h3Index.h b/src/h3lib/include/h3Index.h index a90574e1f8..74773d8507 100644 --- a/src/h3lib/include/h3Index.h +++ b/src/h3lib/include/h3Index.h @@ -21,6 +21,7 @@ #define H3INDEX_H #include "faceijk.h" +#include "geodesicCellBoundary.h" #include "h3api.h" // define's of constants and macros for bitwise manipulation of H3Index's. @@ -181,5 +182,6 @@ DECLSPEC H3Index _zeroIndexDigits(H3Index h, int start, int end); H3Error vec3ToCell(const Vec3d *v, int res, H3Index *out); H3Error cellToVec3(H3Index h3, Vec3d *v); +H3Error cellToGeodesicBoundary(H3Index h3, GeodesicCellBoundary *cb); #endif diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index c88ec98183..9090aea0c0 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -186,6 +186,32 @@ typedef enum { CONTAINMENT_INVALID = 4 ///< This mode is invalid and should not be used } ContainmentMode; +/** + * Bit mask for extracting or setting the containment mode portion of the + * `polygonToCells` and `polygonToCellsExperimental` flags. The four lowest bits + * are reserved for a `ContainmentMode` value. + */ +#define FLAG_CONTAINMENT_MODE_MASK ((uint32_t)(15)) + +#define FLAG_GET_CONTAINMENT_MODE(flags) ((flags)&FLAG_CONTAINMENT_MODE_MASK) + +/** + * Bit mask that enables geodesic polygon filling. Set this bit on the flags to + * have polygon edges follow great-circle arcs instead of planar interpolation. + * + * Geodesic coverage is intended for very large polygons (hundreds of + * kilometers). It is substantially slower than the default planar mode, so + * consider lowering the target resolution when enabling it. + * `CONTAINMENT_CENTER`, `CONTAINMENT_FULL`, and + * `CONTAINMENT_OVERLAPPING` containment modes are valid while the geodesic + * bit is set. + */ +#define FLAG_GEODESIC_MASK ((uint32_t)(16)) + +#define FLAG_GET_GEODESIC(flags) ((flags)&FLAG_GEODESIC_MASK) +#define FLAG_SET_GEODESIC(flags) ((flags) |= FLAG_GEODESIC_MASK) +#define FLAG_RESET_GEODESIC(flags) ((flags) &= ~FLAG_GEODESIC_MASK) + /** @struct LinkedLatLng * @brief A coordinate node in a linked geo structure, part of a linked list */ diff --git a/src/h3lib/include/polyfill.h b/src/h3lib/include/polyfill.h index eccc62e0d2..5e409c52bb 100644 --- a/src/h3lib/include/polyfill.h +++ b/src/h3lib/include/polyfill.h @@ -1,5 +1,5 @@ /* - * Copyright 2023 Uber Technologies, Inc. + * Copyright 2023-2026 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. @@ -22,49 +22,7 @@ #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); +#include "polyfillIterator.h" typedef struct { H3Index cell; // current value diff --git a/src/h3lib/include/polyfillIterator.h b/src/h3lib/include/polyfillIterator.h new file mode 100644 index 0000000000..333c68f455 --- /dev/null +++ b/src/h3lib/include/polyfillIterator.h @@ -0,0 +1,81 @@ +/* + * Copyright 2026 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 polyfillIterator.h + * @brief Iterators used by the polygon-to-cells functionality + */ + +#ifndef POLYFILL_ITERATORS_H +#define POLYFILL_ITERATORS_H + +#include + +#include "bbox.h" +#include "geodesicPolygonInternal.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 responsibility 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 + GeodesicPolygon *geodesicPoly; // Optional geodesic acceleration structure +} IterCellsPolygonCompact; + +DECLSPEC IterCellsPolygonCompact +iterInitPolygonCompact(const GeoPolygon *polygon, int res, uint32_t flags); +DECLSPEC void iterStepPolygonCompact(IterCellsPolygonCompact *iter); +DECLSPEC void iterDestroyPolygonCompact(IterCellsPolygonCompact *iter); +DECLSPEC void iterErrorPolygonCompact(IterCellsPolygonCompact *iter, + H3Error error); + +/** + * Advance to the next cell in the depth-first traversal of the global grid. + * + * @param cell The current cell to advance from. + * @return The successor cell, or `H3_NULL` when the traversal is exhausted. + */ +H3Index nextCell(H3Index cell); + +#endif diff --git a/src/h3lib/include/polygon.h b/src/h3lib/include/polygon.h index 725b206729..2c6d3511ea 100644 --- a/src/h3lib/include/polygon.h +++ b/src/h3lib/include/polygon.h @@ -40,10 +40,6 @@ /** Macro: Whether a GeoLoop is empty */ #define IS_EMPTY_GEOFENCE(geoloop) geoloop->numVerts == 0 -// 1s in the 4 bits defining the polyfill containment mode, 0s elsewhere -#define FLAG_CONTAINMENT_MODE_MASK ((uint32_t)(15)) -#define FLAG_GET_CONTAINMENT_MODE(flags) (flags & FLAG_CONTAINMENT_MODE_MASK) - // Defined in algos.c: void destroyGeoLoop(GeoLoop *loop); void destroyGeoPolygon(GeoPolygon *poly); diff --git a/src/h3lib/include/sphereCapTables.h b/src/h3lib/include/sphereCapTables.h new file mode 100644 index 0000000000..48da1054f9 --- /dev/null +++ b/src/h3lib/include/sphereCapTables.h @@ -0,0 +1,54 @@ +/* + * Copyright 2026 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 sphereCapTables.h + * @brief Shared lookup tables for SphereCap calculations. + * + * These values are used by the library implementation and tests. They are + * defined in a header as `static const` arrays so that including translation + * units receive private copies without exporting additional symbols. + */ + +#ifndef SPHERE_CAP_TABLES_H +#define SPHERE_CAP_TABLES_H + +#include "constants.h" + +/** + * Precomputed cosine values for the bounding cap radius at each resolution. + * Each entry is an IEEE-754 hex float literal equal to + * cos(MAX_EDGE_LENGTH_RADS[res] * CAP_SCALE_FACTOR) + * + * CAP_SCALE_FACTOR (1.01) is the smallest multiplier such that a circle with + * radius MAX_EDGE_LENGTH_RADS[res] * CAP_SCALE_FACTOR fully encloses the cell + * and all of its descendants. + */ +#define H3_SPHERE_CAP_SCALE_FACTOR ((double)1.01) + +static const 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}; + +static const double PRECOMPUTED_COS_RADIUS[MAX_H3_RES + 1] = { + 0x1.f3e3c5bae6cb9p-1, 0x1.fe32bec116ad3p-1, 0x1.ffbdb6ba770aep-1, + 0x1.fff685c0bd115p-1, 0x1.fffea559e01f7p-1, 0x1.ffffce7a554acp-1, + 0x1.fffff8ece6761p-1, 0x1.fffffefd45789p-1, 0x1.ffffffdb09ec8p-1, + 0x1.fffffffab8466p-1, 0x1.ffffffff3ee58p-1, 0x1.ffffffffe469fp-1, + 0x1.fffffffffc0f2p-1, 0x1.ffffffffff6fep-1, 0x1.ffffffffffeb7p-1, + 0x1.fffffffffffd1p-1}; + +#endif diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index aefc701cf2..f1b1811d24 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -877,6 +877,10 @@ H3Error H3_EXPORT(maxPolygonToCellsSize)(const GeoPolygon *geoPolygon, int res, if (flagErr) { return flagErr; } + if (FLAG_GET_GEODESIC(flags)) { + // Legacy planar algorithm does not support geodesic traversal. + return E_OPTION_INVALID; + } // Get the bounding box for the GeoJSON-like struct BBox bbox; const GeoLoop geoloop = geoPolygon->geoloop; @@ -995,6 +999,10 @@ H3Error H3_EXPORT(polygonToCells)(const GeoPolygon *geoPolygon, int res, if (flagErr) { return flagErr; } + if (FLAG_GET_GEODESIC(flags)) { + // Legacy planar algorithm does not support geodesic traversal. + return E_OPTION_INVALID; + } // One of the goals of the polygonToCells algorithm is that two adjacent // polygons with zero overlap have zero overlapping hexagons. That the // hexagons are uniquely assigned. There are a few approaches to take here, diff --git a/src/h3lib/lib/bbox.c b/src/h3lib/lib/bbox.c index 7ea9639297..9bf74cc067 100644 --- a/src/h3lib/lib/bbox.c +++ b/src/h3lib/lib/bbox.c @@ -1,5 +1,5 @@ /* - * Copyright 2016-2021 Uber Technologies, Inc. + * Copyright 2016-2021, 2026 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. @@ -26,6 +26,8 @@ #include "constants.h" #include "h3Index.h" #include "latLng.h" +#include "sphereCapTables.h" +#include "vec3d.h" /** * Width of the bounding box, in rads @@ -309,3 +311,80 @@ void bboxNormalization(const BBox *a, const BBox *b, : aToBTrendsEast ? NORMALIZE_WEST : NORMALIZE_EAST; } + +void aabbEmptyInverted(AABB *box) { + box->min.x = box->min.y = box->min.z = 1.0; + box->max.x = box->max.y = box->max.z = -1.0; +} + +void aabbUpdateWithVec3d(AABB *aabb, const Vec3d *v) { + aabb->min.x = fmin(aabb->min.x, v->x); + aabb->min.y = fmin(aabb->min.y, v->y); + aabb->min.z = fmin(aabb->min.z, v->z); + aabb->max.x = fmax(aabb->max.x, v->x); + aabb->max.y = fmax(aabb->max.y, v->y); + aabb->max.z = fmax(aabb->max.z, v->z); +} + +static inline bool _isOnArcRobust(const Vec3d *p, const Vec3d *v1, + const Vec3d *v2, const Vec3d *n) { + Vec3d c1 = vec3Cross(*v1, *p); + Vec3d c2 = vec3Cross(*p, *v2); + // If p is on the arc, the cross products of (v1, p) and (p, v2) + // should be in the same direction as the normal vector n. + // Using a small tolerance for floating point errors. + return vec3Dot(c1, *n) >= -EPSILON && vec3Dot(c2, *n) >= -EPSILON; +} + +void aabbUpdateWithArcExtrema(AABB *aabb, const Vec3d *v1, const Vec3d *v2, + const Vec3d *n) { + // It is assumed that the AABB has already been updated with the arc's + // endpoints (v1 and v2) before this function is called. + + // Normalize the normal vector for stable projection calculations. + // If the normal is near-zero (degenerate edge, e.g. v1 == v2), there are + // no meaningful extrema to find — the AABB from the endpoints suffices. + Vec3d normal = *n; + if (vec3NormSq(normal) < EPSILON * EPSILON) { + return; + } + vec3Normalize(&normal); + + static const Vec3d UNIT_AXES[3] = { + {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; + + for (int i = 0; i < 3; ++i) { + const Vec3d *axis = &UNIT_AXES[i]; + + // Project the current Cartesian axis onto the great-circle plane. The + // projected direction and its opposite are where that coordinate is + // extremized on the circle. + Vec3d axisCrossNormal = vec3Cross(*axis, normal); + Vec3d extremum_dir = vec3Cross(normal, axisCrossNormal); + double extremum_mag_sq = vec3Dot(extremum_dir, extremum_dir); + + // If the axis is nearly aligned with the normal, there's no unique + // extremum to find. + if (extremum_mag_sq < EPSILON * EPSILON) { + continue; + } + vec3Normalize(&extremum_dir); + + if (_isOnArcRobust(&extremum_dir, v1, v2, &normal)) { + aabbUpdateWithVec3d(aabb, &extremum_dir); + } + + // Check the opposite point on the great circle as well. + Vec3d neg_extremum_dir = {-extremum_dir.x, -extremum_dir.y, + -extremum_dir.z}; + if (_isOnArcRobust(&neg_extremum_dir, v1, v2, &normal)) { + aabbUpdateWithVec3d(aabb, &neg_extremum_dir); + } + } +} + +/** Create a bounding sphere cap for a cell. */ +H3Error cellToSphereCap(H3Index cell, SphereCap *out) { + out->cosRadius = PRECOMPUTED_COS_RADIUS[H3_GET_RESOLUTION(cell)]; + return cellToVec3(cell, &out->center); +} diff --git a/src/h3lib/lib/faceijk.c b/src/h3lib/lib/faceijk.c index 79434d62f2..a87af1a63f 100644 --- a/src/h3lib/lib/faceijk.c +++ b/src/h3lib/lib/faceijk.c @@ -627,6 +627,117 @@ void _faceIjkPentToCellBoundary(const FaceIJK *h, int res, int start, } } +/** + * Generates the cell boundary in 3D coordinates for a pentagonal cell + * given by a FaceIJK address at a specified resolution. + * + * @param h The FaceIJK address of the pentagonal cell. + * @param res The H3 resolution of the cell. + * @param start The first topological vertex to return. + * @param length The number of topological vertexes to return. + * @param g The 3D coordinates of the cell boundary. + */ +void _faceIjkPentToCellBoundaryGeodesic(const FaceIJK *h, int res, int start, + int length, GeodesicCellBoundary *g) { + int adjRes = res; + FaceIJK centerIJK = *h; + FaceIJK fijkVerts[NUM_PENT_VERTS]; + _faceIjkPentToVerts(¢erIJK, &adjRes, fijkVerts); + + // If we're returning the entire loop, we need one more iteration in case + // of a distortion vertex on the last edge + int additionalIteration = length == NUM_PENT_VERTS ? 1 : 0; + + // convert each vertex to lat/lng + // adjust the face of each vertex as appropriate and introduce + // edge-crossing vertices as needed + g->numVerts = 0; + FaceIJK lastFijk = {0}; + for (int vert = start; vert < start + length + additionalIteration; + vert++) { + int v = vert % NUM_PENT_VERTS; + + FaceIJK fijk = fijkVerts[v]; + + _adjustPentVertOverage(&fijk, adjRes); + + // all Class III pentagon edges cross icosa edges + // note that Class II pentagons have vertices on the edge, + // not edge intersections + if (isResolutionClassIII(res) && vert > start) { + // find hex2d of the two vertexes on the last face + + FaceIJK tmpFijk = fijk; + + Vec2d orig2d0; + _ijkToHex2d(&lastFijk.coord, &orig2d0); + + int currentToLastDir = adjacentFaceDir[tmpFijk.face][lastFijk.face]; + + const FaceOrientIJK *fijkOrient = + &faceNeighbors[tmpFijk.face][currentToLastDir]; + + tmpFijk.face = fijkOrient->face; + CoordIJK *ijk = &tmpFijk.coord; + + // rotate and translate for adjacent face + for (int i = 0; i < fijkOrient->ccwRot60; i++) _ijkRotate60ccw(ijk); + + CoordIJK transVec = fijkOrient->translate; + _ijkScale(&transVec, unitScaleByCIIres[adjRes] * 3); + _ijkAdd(ijk, &transVec, ijk); + _ijkNormalize(ijk); + + Vec2d orig2d1; + _ijkToHex2d(ijk, &orig2d1); + + // find the appropriate icosa face edge vertexes + int maxDim = maxDimByCIIres[adjRes]; + Vec2d v0 = {3.0 * maxDim, 0.0}; + Vec2d v1 = {-1.5 * maxDim, 3.0 * M_SQRT3_2 * maxDim}; + Vec2d v2 = {-1.5 * maxDim, -3.0 * M_SQRT3_2 * maxDim}; + + Vec2d *edge0; + Vec2d *edge1; + switch (adjacentFaceDir[tmpFijk.face][fijk.face]) { + case IJ: + edge0 = &v0; + edge1 = &v1; + break; + case JK: + edge0 = &v1; + edge1 = &v2; + break; + case KI: + default: + assert(adjacentFaceDir[tmpFijk.face][fijk.face] == KI); + edge0 = &v2; + edge1 = &v0; + break; + } + + // find the intersection and add the lat/lng point to the result + Vec2d inter; + _v2dIntersect(&orig2d0, &orig2d1, edge0, edge1, &inter); + _hex2dToVec3(&inter, tmpFijk.face, adjRes, 1, + &g->verts[g->numVerts]); + g->numVerts++; + } + + // convert vertex to lat/lng and add to the result + // vert == start + NUM_PENT_VERTS is only used to test for possible + // intersection on last edge + if (vert < start + NUM_PENT_VERTS) { + Vec2d vec; + _ijkToHex2d(&fijk.coord, &vec); + _hex2dToVec3(&vec, fijk.face, adjRes, 1, &g->verts[g->numVerts]); + g->numVerts++; + } + + lastFijk = fijk; + } +} + /** * Get the vertices of a pentagon cell as substrate FaceIJK addresses * @@ -804,6 +915,119 @@ void _faceIjkToCellBoundary(const FaceIJK *h, int res, int start, int length, } } +/** + * Generates the cell boundary in 3D coordinates for a cell given by a + * FaceIJK address at a specified resolution. + * + * @param h The FaceIJK address of the cell. + * @param res The H3 resolution of the cell. + * @param start The first topological vertex to return. + * @param length The number of topological vertexes to return. + * @param g The 3D coordinates of the cell boundary. + */ +void _faceIjkToCellBoundaryGeodesic(const FaceIJK *h, int res, int start, + int length, GeodesicCellBoundary *g) { + int adjRes = res; + FaceIJK centerIJK = *h; + FaceIJK fijkVerts[NUM_HEX_VERTS]; + _faceIjkToVerts(¢erIJK, &adjRes, fijkVerts); + + // If we're returning the entire loop, we need one more iteration in case + // of a distortion vertex on the last edge + int additionalIteration = length == NUM_HEX_VERTS ? 1 : 0; + + // convert each vertex to xyz + // adjust the face of each vertex as appropriate and introduce + // edge-crossing vertices as needed + g->numVerts = 0; + int lastFace = -1; + Overage lastOverage = NO_OVERAGE; + for (int vert = start; vert < start + length + additionalIteration; + vert++) { + int v = vert % NUM_HEX_VERTS; + + FaceIJK fijk = fijkVerts[v]; + + const int pentLeading4 = 0; + Overage overage = _adjustOverageClassII(&fijk, adjRes, pentLeading4, 1); + + /* + Check for edge-crossing. Each face of the underlying icosahedron is a + different projection plane. So if an edge of the hexagon crosses an + icosahedron edge, an additional vertex must be introduced at that + intersection point. Then each half of the cell edge can be projected + to geographic coordinates using the appropriate icosahedron face + projection. Note that Class II cell edges have vertices on the face + edge, with no edge line intersections. + */ + if (isResolutionClassIII(res) && vert > start && + fijk.face != lastFace && lastOverage != FACE_EDGE) { + // find hex2d of the two vertexes on original face + int lastV = (v + 5) % NUM_HEX_VERTS; + Vec2d orig2d0; + _ijkToHex2d(&fijkVerts[lastV].coord, &orig2d0); + + Vec2d orig2d1; + _ijkToHex2d(&fijkVerts[v].coord, &orig2d1); + + // find the appropriate icosa face edge vertexes + int maxDim = maxDimByCIIres[adjRes]; + Vec2d v0 = {3.0 * maxDim, 0.0}; + Vec2d v1 = {-1.5 * maxDim, 3.0 * M_SQRT3_2 * maxDim}; + Vec2d v2 = {-1.5 * maxDim, -3.0 * M_SQRT3_2 * maxDim}; + + int face2 = ((lastFace == centerIJK.face) ? fijk.face : lastFace); + Vec2d *edge0; + Vec2d *edge1; + switch (adjacentFaceDir[centerIJK.face][face2]) { + case IJ: + edge0 = &v0; + edge1 = &v1; + break; + case JK: + edge0 = &v1; + edge1 = &v2; + break; + // case KI: + default: + assert(adjacentFaceDir[centerIJK.face][face2] == KI); + edge0 = &v2; + edge1 = &v0; + break; + } + + // find the intersection and add the lat/lng point to the result + Vec2d inter; + _v2dIntersect(&orig2d0, &orig2d1, edge0, edge1, &inter); + /* + If a point of intersection occurs at a hexagon vertex, then each + adjacent hexagon edge will lie completely on a single icosahedron + face, and no additional vertex is required. + */ + bool isIntersectionAtVertex = _v2dAlmostEquals(&orig2d0, &inter) || + _v2dAlmostEquals(&orig2d1, &inter); + if (!isIntersectionAtVertex) { + _hex2dToVec3(&inter, centerIJK.face, adjRes, 1, + &g->verts[g->numVerts]); + g->numVerts++; + } + } + + // convert vertex to lat/lng and add to the result + // vert == start + NUM_HEX_VERTS is only used to test for possible + // intersection on last edge + if (vert < start + NUM_HEX_VERTS) { + Vec2d vec; + _ijkToHex2d(&fijk.coord, &vec); + _hex2dToVec3(&vec, fijk.face, adjRes, 1, &g->verts[g->numVerts]); + g->numVerts++; + } + + lastFace = fijk.face; + lastOverage = overage; + } +} + /** * Get the vertices of a cell as substrate FaceIJK addresses * diff --git a/src/h3lib/lib/geodesicIterator.c b/src/h3lib/lib/geodesicIterator.c new file mode 100644 index 0000000000..6c3b51c2b2 --- /dev/null +++ b/src/h3lib/lib/geodesicIterator.c @@ -0,0 +1,198 @@ +/* + * Copyright 2026 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 geodesicIterator.c + * @brief Manages geodesic polygon iterator lifecycle and traversal helpers. + */ + +#include "geodesicIterator.h" + +#include "geodesicPolygonInternal.h" +#include "h3Index.h" + +#define H3_CHECK(h3_call, iter_ptr) \ + do { \ + H3Error _err = (h3_call); \ + if (_err != E_SUCCESS) { \ + iterErrorPolygonCompact((iter_ptr), _err); \ + return; \ + } \ + } while (0) + +static GeodesicPolygon *_getOrCreateGeodesicPolygon( + IterCellsPolygonCompact *iter) { + if (iter->geodesicPoly == NULL) { + if (iter->_polygon->geoloop.numVerts < 3) { + iterErrorPolygonCompact(iter, E_DOMAIN); + return NULL; + } + if (iter->_polygon->numHoles > 0 && iter->_polygon->holes == NULL) { + iterErrorPolygonCompact(iter, E_DOMAIN); + return NULL; + } + for (int i = 0; i < iter->_polygon->numHoles; i++) { + if (iter->_polygon->holes[i].numVerts < 3) { + iterErrorPolygonCompact(iter, E_DOMAIN); + return NULL; + } + } + GeodesicPolygon *poly = NULL; + H3Error createErr = geodesicPolygonCreate(iter->_polygon, &poly); + if (createErr != E_SUCCESS) { + iterErrorPolygonCompact(iter, createErr); + return NULL; + } + iter->geodesicPoly = poly; + } + return iter->geodesicPoly; +} + +void geodesicIteratorDestroyState(IterCellsPolygonCompact *iter) { + if (!iter || !iter->geodesicPoly) { + return; + } + geodesicPolygonDestroy(iter->geodesicPoly); + iter->geodesicPoly = NULL; +} + +void geodesicIteratorStep(IterCellsPolygonCompact *iter, H3Index cell) { + GeodesicPolygon *poly = _getOrCreateGeodesicPolygon(iter); + if (!poly) { + return; + } + + ContainmentMode mode = FLAG_GET_CONTAINMENT_MODE(iter->_flags); + + while (cell) { + int cellRes = H3_GET_RESOLUTION(cell); + + // Coarse cells: prune by cap/AABB, otherwise descend into children. + SphereCap cap; + H3_CHECK(cellToSphereCap(cell, &cap), iter); + if (cellRes < iter->_res) { + if (!geodesicPolygonCapIntersects(poly, &cap)) { + cell = nextCell(cell); + continue; + } + + H3Index child; + H3_CHECK(H3_EXPORT(cellToCenterChild)(cell, cellRes + 1, &child), + iter); + cell = child; + continue; + } + + // Fine cells: start with center-in-polygon classification. + Vec3d cellCenter; + H3_CHECK(cellToVec3(cell, &cellCenter), iter); + const bool pointInside = + geodesicPolygonContainsPoint(poly, &cellCenter); + + // CENTER mode uses center-only semantics. + if (mode == CONTAINMENT_CENTER) { + if (pointInside) { + iter->cell = cell; + return; + } + cell = nextCell(cell); + continue; + } + + // OVERLAPPING can accept immediately when center is inside. + if (mode == CONTAINMENT_OVERLAPPING && pointInside) { + iter->cell = cell; + return; + } + + // FULL cannot match when the center is outside. + if (mode == CONTAINMENT_FULL && !pointInside) { + cell = nextCell(cell); + continue; + } + + // Expensive boundary tests are only needed for unresolved cases. + GeodesicCellBoundary boundary; + H3_CHECK(cellToGeodesicBoundary(cell, &boundary), iter); + + const bool boundaryIntersection = + geodesicPolygonBoundaryIntersects(poly, &boundary, &cap); + + if (mode == CONTAINMENT_FULL) { + // A cell is fully contained only when all of the following hold + // 1. the cell boundary does not cross the polygon boundary + // (outer loop or any hole), + // 2. a cell boundary vertex lies inside the filled polygon, and + // 3. no hole lies entirely inside the cell. + if (boundaryIntersection) { + cell = nextCell(cell); + continue; + } + if (!geodesicPolygonContainsPoint(poly, &boundary.verts[0])) { + cell = nextCell(cell); + continue; + } + // A hole whose first vertex lands in this cell lies inside it + bool holeInsideCell = false; + for (int hi = 0; hi < poly->numHoles; hi++) { + H3Index holeCell; + H3_CHECK(vec3ToCell(&poly->holes[hi].edges[0].vert, cellRes, + &holeCell), + iter); + if (holeCell == cell) { + holeInsideCell = true; + break; + } + } + if (holeInsideCell) { + cell = nextCell(cell); + continue; + } + iter->cell = cell; + return; + } + + // OVERLAPPING and center outside - match if we intersect + if (boundaryIntersection) { + iter->cell = cell; + return; + } + + // Center is outside and no edge intersection. The polygon may + // be inside the cell. Check if any loop vertex (outer shell or hole) + // falls inside the cell. An outer vertex here means the polygon is + // smaller than the cell; a hole vertex here means the hole is inside + // the cell, so the filled region (outer minus holes) still overlaps it. + H3Index polygonCell; + Vec3d *firstVert = &poly->geoloop.edges[0].vert; + H3_CHECK(vec3ToCell(firstVert, cellRes, &polygonCell), iter); + if (polygonCell == cell) { + iter->cell = cell; + return; + } + for (int hi = 0; hi < poly->numHoles; hi++) { + Vec3d *holeVert = &poly->holes[hi].edges[0].vert; + H3_CHECK(vec3ToCell(holeVert, cellRes, &polygonCell), iter); + if (polygonCell == cell) { + iter->cell = cell; + return; + } + } + + cell = nextCell(cell); + } + + iterDestroyPolygonCompact(iter); +} diff --git a/src/h3lib/lib/geodesicPolygon.c b/src/h3lib/lib/geodesicPolygon.c new file mode 100644 index 0000000000..45a94d6008 --- /dev/null +++ b/src/h3lib/lib/geodesicPolygon.c @@ -0,0 +1,576 @@ +/* + * Copyright 2026 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 geodesicPolygon.c + * @brief Implements geodesic polygon conversions and spatial predicates. + */ + +#include + +#include "alloc.h" +#include "constants.h" + +// Tolerance for classifying arc endpoints as being on the same side of a +// great-circle plane. See _geodesicEdgesCross for derivation. +static const double SIDE_EPS = 1e-28; +#include "geodesicPolygonInternal.h" +#include "h3api.h" + +/** + * Determine whether two geodesic edges intersect on the unit sphere. + * + * The test classifies edge endpoints against each other's great-circle planes + * to cull obvious misses, handles near-parallel arcs with a 1-D projection + * fallback, and otherwise checks that the great-circle intersection lies + * between both segments. + */ +static bool _geodesicEdgesCross(const Vec3d *a1, const Vec3d *a2, + const Vec3d *b1, const Vec3d *b2, + const Vec3d *normalB) { + Vec3d normalA; + normalA = vec3Cross(*a1, *a2); + + const double b1Side = vec3Dot(normalA, *b1); + const double b2Side = vec3Dot(normalA, *b2); + const double a1Side = vec3Dot(*normalB, *a1); + const double a2Side = vec3Dot(*normalB, *a2); + + // For coincident arcs (shared cell/polygon boundary edges), the triple + // scalar products b1Side/b2Side are mathematically zero but may reach + // ~2e-32 in FP due to per-operand rounding (~1e-16 per dot product, two + // of them multiplied together). EPSILON*EPSILON = 1e-32 sits right at + // that noise floor, so some coincident pairs exceed it and are wrongly + // rejected. SIDE_EPS sits four orders of magnitude above + // machine_epsilon^2 (~5e-32) but far below the smallest genuine same-side + // product observed in practice (~1e-18) — so coincident arcs reliably + // fall through to the collinear-overlap check while non-coincident + // same-side pairs still early-exit. + if ((b1Side * b2Side > SIDE_EPS) || (a1Side * a2Side > SIDE_EPS)) { + return false; + } + + Vec3d intersectionLine; + intersectionLine = vec3Cross(normalA, *normalB); + + if (vec3NormSq(intersectionLine) < EPSILON * EPSILON) { + const Vec3d refDir = {a2->x - a1->x, a2->y - a1->y, a2->z - a1->z}; + if (vec3NormSq(refDir) < EPSILON * EPSILON) { + return false; + } + + double projA1 = vec3Dot(*a1, refDir); + double projA2 = vec3Dot(*a2, refDir); + double projB1 = vec3Dot(*b1, refDir); + double projB2 = vec3Dot(*b2, refDir); + + if (projA1 > projA2) { + double tmp = projA1; + projA1 = projA2; + projA2 = tmp; + } + if (projB1 > projB2) { + double tmp = projB1; + projB1 = projB2; + projB2 = tmp; + } + + return (projA1 <= projB2) && (projB1 <= projA2); + } + + const Vec3d midASum = {a1->x + a2->x, a1->y + a2->y, a1->z + a2->z}; + const Vec3d midBSum = {b1->x + b2->x, b1->y + b2->y, b1->z + b2->z}; + + if (vec3NormSq(midASum) < EPSILON * EPSILON || + vec3NormSq(midBSum) < EPSILON * EPSILON) { + return true; + } + + double dotA = vec3Dot(intersectionLine, midASum); + double dotB = vec3Dot(intersectionLine, midBSum); + + return (dotA * dotB >= -EPSILON); +} + +/** + * Test whether a point lies inside a geodesic loop using spherical winding. + * + * The algorithm walks the loop once, accumulating the signed angle subtended by + * consecutive vertices as seen from the query point. When the total winding + * exceeds pi in magnitude the point is inside; otherwise it is outside. + * + * A hemisphere-pole early rejection ensures that points on the far side of the + * sphere (where the winding formula can give spurious +-2pi results) are + * correctly classified as outside. This is safe because the polygon must fit + * within a hemisphere; any point in the opposite hemisphere is guaranteed to be + * outside. + */ +static bool _geodesicLoopContainsPoint(const GeodesicLoop *loop, + const Vec3d *pointVec) { + // Early rejection: if the loop has a known hemisphere pole and the point + // lies in the opposite hemisphere, it cannot be inside the polygon. + // This also prevents spurious winding-number results for antipodal points. + if (loop->hasHemispherePole && + vec3Dot(loop->hemispherePole, *pointVec) < -1e-10) { + return false; + } + + double totalAngle = 0; + double dotPV1 = vec3Dot(*pointVec, loop->edges[0].vert); + + for (int i = 0; i < loop->numEdges; i++) { + Vec3d *vert2 = &loop->edges[(i + 1) % loop->numEdges].vert; + double dotPV2 = vec3Dot(*pointVec, *vert2); + + double y = vec3Dot(*pointVec, loop->edges[i].edgeCross); + double x = loop->edges[i].edgeDot - dotPV1 * dotPV2; + + totalAngle += atan2(y, x); + + dotPV1 = dotPV2; + } + + return fabs(totalAngle) > M_PI; +} + +/** + * Quickly reject a sphere cap and AABB that cannot intersect. + * Returns false ONLY when intersection is definitively impossible. + */ +static bool _geodesicSphereCapOverlapsAABB(const SphereCap *cap, + const AABB *aabb) { + // Cos comparisons require more accuracy - use bigger epsilon + const double epsilon = 1e-8; + + // 1. Far point test - checks if the farthest AABB corner + // falls outside the cap + Vec3d farPoint; + farPoint.x = cap->center.x >= 0 ? aabb->max.x : aabb->min.x; + farPoint.y = cap->center.y >= 0 ? aabb->max.y : aabb->min.y; + farPoint.z = cap->center.z >= 0 ? aabb->max.z : aabb->min.z; + + double farDot = vec3Dot(farPoint, cap->center); + if (farDot < cap->cosRadius - epsilon) { + return false; + } + + // 2. Near origin test - checks if the closest AABB point + // is outside the unit sphere + Vec3d nearOrigin; + nearOrigin.x = fmax(aabb->min.x, fmin(aabb->max.x, 0.0)); + nearOrigin.y = fmax(aabb->min.y, fmin(aabb->max.y, 0.0)); + nearOrigin.z = fmax(aabb->min.z, fmin(aabb->max.z, 0.0)); + + double distSq = vec3Dot(nearOrigin, nearOrigin); + if (distSq > 1.0 + epsilon) { + return false; + } + + return true; +} + +static void _geodesicLoopToAABB(const GeodesicLoop *loop, AABB *out) { + aabbEmptyInverted(out); + + for (int i = 0; i < loop->numEdges; i++) { + aabbUpdateWithVec3d(out, &loop->edges[i].vert); + + aabbUpdateWithArcExtrema(out, &loop->edges[i].vert, + &loop->edges[(i + 1) % loop->numEdges].vert, + &loop->edges[i].edgeCross); + } + + // Probe cardinal axes and expand the box if needed + Vec3d testVecs[] = {{1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, + {0, -1, 0}, {0, 0, 1}, {0, 0, -1}}; + for (int i = 0; i < 6; i++) { + if (_geodesicLoopContainsPoint(loop, &testVecs[i])) { + aabbUpdateWithVec3d(out, &testVecs[i]); + } + } +} + +static bool _candidateContainsVerts(const Vec3d *candidate, + const GeodesicEdge *edges, int numEdges) { + if (vec3NormSq(*candidate) < EPSILON * EPSILON) { + return false; + } + + Vec3d normalized = *candidate; + vec3Normalize(&normalized); + + const double hemisphereEpsilon = -1e-12; + for (int i = 0; i < numEdges; i++) { + if (vec3Dot(normalized, edges[i].vert) < hemisphereEpsilon) { + return false; + } + } + return true; +} + +/** + * Test whether all loop vertices fit within some closed hemisphere. + * + * A hemisphere can be described by a pole vector c such that dot(c, v) >= 0 + * for every vertex v in the loop. + * + * This checks a deterministic set of candidate poles and returns true once + * one contains all vertices: + * - The vertex centroid + * - Each vertex + * - cross(v_i, v_j) and its negation for each vertex pair + * + * In the common case, the centroid check succeeds in O(n). If not, the + * fallback pairwise search is more expensive (O(n^3) worst case). + */ +static bool _geoLoopVerticesFitHemisphere(const GeodesicEdge *edges, + int numEdges, Vec3d *poleOut) { + Vec3d centroid = {0}; + for (int i = 0; i < numEdges; i++) { + centroid.x += edges[i].vert.x; + centroid.y += edges[i].vert.y; + centroid.z += edges[i].vert.z; + } + + // Near-zero centroid direction is numerically inconclusive, so skip it. + // More tight epsilon check for more conclusive result + if (vec3NormSq(centroid) >= 1e-12) { + if (_candidateContainsVerts(¢roid, edges, numEdges)) { + if (poleOut) { + *poleOut = centroid; + vec3Normalize(poleOut); + } + return true; + } + } + + for (int i = 0; i < numEdges; i++) { + if (_candidateContainsVerts(&edges[i].vert, edges, numEdges)) { + if (poleOut) { + *poleOut = edges[i].vert; + } + return true; + } + } + + for (int i = 0; i < numEdges; i++) { + for (int j = i + 1; j < numEdges; j++) { + Vec3d candidate; + candidate = vec3Cross(edges[i].vert, edges[j].vert); + + if (_candidateContainsVerts(&candidate, edges, numEdges)) { + if (poleOut) { + *poleOut = candidate; + vec3Normalize(poleOut); + } + return true; + } + + Vec3d opposite = { + .x = -candidate.x, .y = -candidate.y, .z = -candidate.z}; + if (_candidateContainsVerts(&opposite, edges, numEdges)) { + if (poleOut) { + *poleOut = opposite; + vec3Normalize(poleOut); + } + return true; + } + } + } + + return false; +} + +/** + * Test whether all loop vertices lie on a single great circle. + */ +static bool _geoLoopVerticesOnGreatCircle(const GeodesicEdge *edges, + int numEdges) { + // Keep legacy handling for point/line/triangle degeneracies and only + // reject larger boundary-only loops that produce ambiguous containment. + if (numEdges < 4) { + return false; + } + + Vec3d normal = {0}; + bool foundNormal = false; + for (int i = 0; i < numEdges && !foundNormal; i++) { + for (int j = i + 1; j < numEdges; j++) { + normal = vec3Cross(edges[i].vert, edges[j].vert); + if (vec3NormSq(normal) > 1e-12) { + foundNormal = true; + break; + } + } + } + + if (!foundNormal) { + return false; + } + + vec3Normalize(&normal); + const double greatCircleEpsilon = 1e-14; + for (int i = 0; i < numEdges; i++) { + if (fabs(vec3Dot(normal, edges[i].vert)) > greatCircleEpsilon) { + return false; + } + } + + return true; +} + +static H3Error _geodesicLoopFromGeo(const GeoLoop *loop, GeodesicLoop *out, + bool rejectLargeLoop) { + if (!out || !loop || loop->numVerts < 3) { + return E_DOMAIN; + } + *out = (GeodesicLoop){0}; + + const int n = loop->numVerts; + GeodesicEdge *edges = H3_MEMORY(calloc)(n, sizeof(GeodesicEdge)); + if (!edges) { + return E_MEMORY_ALLOC; + } + + out->edges = edges; + out->numEdges = n; + out->centroid.x = 0; + out->centroid.y = 0; + out->centroid.z = 0; + + for (int i = 0; i < n; i++) { + if (!isfinite(loop->verts[i].lat) || !isfinite(loop->verts[i].lng)) { + H3_MEMORY(free)(edges); + out->edges = NULL; + out->numEdges = 0; + return E_DOMAIN; + } + edges[i].vert = latLngToVec3(loop->verts[i]); + } + + // Great-circle boundary loops are always ambiguous (zero enclosed area) + // regardless of polygon size, so reject them unconditionally. + if (_geoLoopVerticesOnGreatCircle(edges, n)) { + H3_MEMORY(free)(edges); + out->edges = NULL; + out->numEdges = 0; + return E_DOMAIN; + } + + if (rejectLargeLoop) { + Vec3d pole; + if (!_geoLoopVerticesFitHemisphere(edges, n, &pole)) { + H3_MEMORY(free)(edges); + out->edges = NULL; + out->numEdges = 0; + return E_DOMAIN; + } + out->hemispherePole = pole; + out->hasHemispherePole = true; + } + + for (int i = 0; i < n; i++) { + const Vec3d *v1 = &edges[i].vert; + const Vec3d *v2 = &edges[(i + 1) % n].vert; + + out->centroid.x += v1->x; + out->centroid.y += v1->y; + out->centroid.z += v1->z; + + edges[i].edgeCross = vec3Cross(*v1, *v2); + edges[i].edgeDot = vec3Dot(*v1, *v2); + + AABB *box = &edges[i].aabb; + aabbEmptyInverted(box); + aabbUpdateWithVec3d(box, v1); + aabbUpdateWithVec3d(box, v2); + aabbUpdateWithArcExtrema(box, v1, v2, &edges[i].edgeCross); + } + + // Normalize the centroid to a unit vector for use in AABB construction + // and other geometric queries. + vec3Normalize(&out->centroid); + + return E_SUCCESS; +} + +static void _geodesicLoopDestroy(GeodesicLoop *loop) { + if (!loop || !loop->edges) { + return; + } + + H3_MEMORY(free)(loop->edges); + loop->edges = NULL; + loop->numEdges = 0; +} + +H3Error geodesicPolygonCreate(const GeoPolygon *polygon, + GeodesicPolygon **out) { + if (!polygon || !out || polygon->geoloop.numVerts < 3) { + return E_DOMAIN; + } + + GeodesicPolygon *result = H3_MEMORY(calloc)(1, sizeof(GeodesicPolygon)); + if (!result) { + return E_MEMORY_ALLOC; + } + + H3Error loopErr = + _geodesicLoopFromGeo(&polygon->geoloop, &result->geoloop, true); + if (loopErr != E_SUCCESS) { + _geodesicLoopDestroy(&result->geoloop); + H3_MEMORY(free)(result); + return loopErr; + } + + const int holeCount = polygon->numHoles; + result->numHoles = holeCount; + if (holeCount > 0) { + if (!polygon->holes) { + _geodesicLoopDestroy(&result->geoloop); + H3_MEMORY(free)(result); + return E_DOMAIN; + } + result->holes = H3_MEMORY(calloc)(holeCount, sizeof(GeodesicLoop)); + if (!result->holes) { + _geodesicLoopDestroy(&result->geoloop); + H3_MEMORY(free)(result); + return E_MEMORY_ALLOC; + } + for (int i = 0; i < holeCount; i++) { + if (polygon->holes[i].numVerts < 3) { + for (int j = 0; j < i; j++) { + _geodesicLoopDestroy(&result->holes[j]); + } + H3_MEMORY(free)(result->holes); + _geodesicLoopDestroy(&result->geoloop); + H3_MEMORY(free)(result); + return E_DOMAIN; + } + H3Error holeErr = _geodesicLoopFromGeo(&polygon->holes[i], + &result->holes[i], false); + if (holeErr != E_SUCCESS) { + for (int j = 0; j < i; j++) { + _geodesicLoopDestroy(&result->holes[j]); + } + H3_MEMORY(free)(result->holes); + _geodesicLoopDestroy(&result->geoloop); + H3_MEMORY(free)(result); + return holeErr; + } + } + } + + _geodesicLoopToAABB(&result->geoloop, &result->aabb); + + *out = result; + return E_SUCCESS; +} + +void geodesicPolygonDestroy(GeodesicPolygon *polygon) { + if (!polygon) { + return; + } + + _geodesicLoopDestroy(&polygon->geoloop); + + if (polygon->holes) { + for (int i = 0; i < polygon->numHoles; i++) { + _geodesicLoopDestroy(&polygon->holes[i]); + } + H3_MEMORY(free)(polygon->holes); + } + + H3_MEMORY(free)(polygon); +} + +bool geodesicPolygonCapIntersects(const GeodesicPolygon *polygon, + const SphereCap *cap) { + if (!polygon || !cap) { + return false; + } + + return _geodesicSphereCapOverlapsAABB(cap, &polygon->aabb); +} + +/** + * Check whether a geodesic polygon boundary intersects a cell boundary. + * + * For each polygon loop (outer shell and holes) we prune arcs whose bounding + * boxes do not overlap the query cap. Remaining arcs are tested pairwise + * against the cell boundary segments with `_geodesicEdgesCross`. Early returns + * avoid unnecessary great-circle computations. + */ +bool geodesicPolygonBoundaryIntersects(const GeodesicPolygon *polygon, + const GeodesicCellBoundary *boundary, + const SphereCap *cap) { + if (!polygon || !boundary || !cap) { + return false; + } + + const GeodesicLoop *loops = &polygon->geoloop; + const int loopCount = 1 + polygon->numHoles; + + for (int loopIdx = 0; loopIdx < loopCount; loopIdx++) { + const GeodesicLoop *loop = + (loopIdx == 0) ? loops : &polygon->holes[loopIdx - 1]; + + for (int i = 0; i < loop->numEdges; i++) { + if (!_geodesicSphereCapOverlapsAABB(cap, &loop->edges[i].aabb)) { + continue; + } + + int nextI = (i + 1) % loop->numEdges; + for (int j = 0; j < boundary->numVerts; j++) { + int nextJ = (j + 1) % boundary->numVerts; + if (_geodesicEdgesCross( + &boundary->verts[j], &boundary->verts[nextJ], + &loop->edges[i].vert, &loop->edges[nextI].vert, + &loop->edges[i].edgeCross)) { + return true; + } + } + } + } + + return false; +} + +/** + * Evaluate whether a 3D unit vector lies inside a polygon with optional holes. + * + * The outer loop is tested for containment first; if it fails we can reject the + * point immediately. When the point is inside the shell we ensure it does not + * land inside any hole loop, mirroring the typical point-in-polygon predicate + * on the sphere. + */ +bool geodesicPolygonContainsPoint(const GeodesicPolygon *polygon, + const Vec3d *point) { + if (!polygon || !point) { + return false; + } + + if (!_geodesicLoopContainsPoint(&polygon->geoloop, point)) { + return false; + } + + for (int i = 0; i < polygon->numHoles; i++) { + if (_geodesicLoopContainsPoint(&polygon->holes[i], point)) { + return false; + } + } + + return true; +} diff --git a/src/h3lib/lib/h3Index.c b/src/h3lib/lib/h3Index.c index 64c11c9ef5..e11918100c 100644 --- a/src/h3lib/lib/h3Index.c +++ b/src/h3lib/lib/h3Index.c @@ -1221,6 +1221,28 @@ H3Error H3_EXPORT(cellToBoundary)(H3Index h3, CellBoundary *cb) { return E_SUCCESS; } +/** + * Determines the cell boundary in spherical coordinates for an H3 index. + * + * @param h3 The H3 index. + * @param cb The boundary of the H3 cell in spherical coordinates. + */ +H3Error cellToGeodesicBoundary(H3Index h3, GeodesicCellBoundary *cb) { + FaceIJK fijk; + H3Error e = _h3ToFaceIjk(h3, &fijk); + if (e) { + return e; + } + if (H3_EXPORT(isPentagon)(h3)) { + _faceIjkPentToCellBoundaryGeodesic(&fijk, H3_GET_RESOLUTION(h3), 0, + NUM_PENT_VERTS, cb); + } else { + _faceIjkToCellBoundaryGeodesic(&fijk, H3_GET_RESOLUTION(h3), 0, + NUM_HEX_VERTS, cb); + } + return E_SUCCESS; +} + /** * Returns the max number of possible icosahedron faces an H3 index * may intersect. diff --git a/src/h3lib/lib/polyfill.c b/src/h3lib/lib/polyfill.c index 40e67ecdcf..a2439cb0aa 100644 --- a/src/h3lib/lib/polyfill.c +++ b/src/h3lib/lib/polyfill.c @@ -1,5 +1,5 @@ /* - * Copyright 2023 Uber Technologies, Inc. + * Copyright 2023-2026 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. @@ -25,9 +25,11 @@ #include "alloc.h" #include "baseCells.h" #include "coordijk.h" +#include "geodesicIterator.h" #include "h3Assert.h" #include "h3Index.h" #include "polygon.h" +#include "sphereCapTables.h" // Factor by which to scale the cell bounding box to include all cells. // This was determined empirically by finding the smallest factor that @@ -39,17 +41,6 @@ // 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, @@ -277,8 +268,7 @@ H3Index baseCellNumToCell(int baseCellNum) { return baseCell; } -static void iterErrorPolygonCompact(IterCellsPolygonCompact *iter, - H3Error error) { +void iterErrorPolygonCompact(IterCellsPolygonCompact *iter, H3Error error) { iterDestroyPolygonCompact(iter); iter->error = error; } @@ -287,7 +277,7 @@ static void iterErrorPolygonCompact(IterCellsPolygonCompact *iter, * Given a cell, find the next cell in the sequence of all cells * to check in the iteration. */ -static H3Index nextCell(H3Index cell) { +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) @@ -332,7 +322,8 @@ static IterCellsPolygonCompact _iterInitPolygonCompact( ._res = res, ._flags = flags, ._bboxes = NULL, - ._started = false}; + ._started = false, + .geodesicPoly = NULL}; if (res < 0 || res > MAX_H3_RES) { iterErrorPolygonCompact(&iter, E_RES_DOMAIN); @@ -345,6 +336,11 @@ static IterCellsPolygonCompact _iterInitPolygonCompact( return iter; } + if (polygon->numHoles > 0 && polygon->holes == NULL) { + iterErrorPolygonCompact(&iter, E_DOMAIN); + 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)); @@ -426,6 +422,11 @@ void iterStepPolygonCompact(IterCellsPolygonCompact *iter) { return; } + if (FLAG_GET_GEODESIC(iter->_flags)) { + geodesicIteratorStep(iter, cell); + return; + } + ContainmentMode mode = FLAG_GET_CONTAINMENT_MODE(iter->_flags); while (cell) { @@ -597,13 +598,14 @@ void iterStepPolygonCompact(IterCellsPolygonCompact *iter) { void iterDestroyPolygonCompact(IterCellsPolygonCompact *iter) { if (iter->_bboxes) { H3_MEMORY(free)(iter->_bboxes); + iter->_bboxes = NULL; } + geodesicIteratorDestroyState(iter); iter->cell = H3_NULL; iter->error = E_SUCCESS; iter->_polygon = NULL; iter->_res = -1; iter->_flags = 0; - iter->_bboxes = NULL; } /** @@ -746,8 +748,16 @@ H3Error H3_EXPORT(maxPolygonToCellsSizeExperimental)(const GeoPolygon *polygon, return iter.error; } - // Ignore the requested flags and use the faster overlapping-bbox mode - iter._flags = CONTAINMENT_OVERLAPPING_BBOX; + // Size estimates must preserve geodesic traversal semantics; otherwise the + // estimate can undercount and violate the API contract. + if (FLAG_GET_GEODESIC(iter._flags)) { + iter._flags = CONTAINMENT_OVERLAPPING; + FLAG_SET_GEODESIC(iter._flags); + } else { + // In planar mode we can safely use overlapping-bbox to speed up + // estimation while remaining an upper bound. + iter._flags = CONTAINMENT_OVERLAPPING_BBOX; + } // Get a (very) rough area of the polygon bounding box BBox *polygonBBox = &iter._bboxes[0]; @@ -756,9 +766,9 @@ H3Error H3_EXPORT(maxPolygonToCellsSizeExperimental)(const GeoPolygon *polygon, cos(fmin(fabs(polygonBBox->north), fabs(polygonBBox->south))) * EARTH_RADIUS_KM * EARTH_RADIUS_KM; - // Determine the res for the size estimate, based on a (very) rough estimate - // of the number of cells at various resolutions that would fit in the - // polygon. All we need here is a general order of magnitude. + // Determine the res for the size estimate, based on a (very) rough + // estimate of the number of cells at various resolutions that would fit in + // the polygon. All we need here is a general order of magnitude. while (iter._res > 0 && polygonBBoxAreaKm2 / getAverageCellArea(iter._res - 1) > MAX_SIZE_CELL_THRESHOLD) { diff --git a/src/h3lib/lib/polygon.c b/src/h3lib/lib/polygon.c index 18f5b858cf..7537617d55 100644 --- a/src/h3lib/lib/polygon.c +++ b/src/h3lib/lib/polygon.c @@ -49,10 +49,31 @@ * @return Whether the flags are valid */ H3Error validatePolygonFlags(uint32_t flags) { - if (flags & (~FLAG_CONTAINMENT_MODE_MASK) || - FLAG_GET_CONTAINMENT_MODE(flags) >= CONTAINMENT_INVALID) { + ContainmentMode containmentMode = FLAG_GET_CONTAINMENT_MODE(flags); + bool isGeodesic = FLAG_GET_GEODESIC(flags); + + // Check containment mode validity based on geodesic flag + if (isGeodesic) { + // For geodesic: CENTER, FULL, and OVERLAPPING are valid + if (containmentMode != CONTAINMENT_CENTER && + containmentMode != CONTAINMENT_FULL && + containmentMode != CONTAINMENT_OVERLAPPING) { + return E_OPTION_INVALID; + } + } else { + // For non-geodesic: everything < CONTAINMENT_INVALID is valid + if (containmentMode >= CONTAINMENT_INVALID) { + return E_OPTION_INVALID; + } + } + + // Check if any invalid flags are set (only containment mode and geodesic + // are valid) + uint32_t validFlags = FLAG_CONTAINMENT_MODE_MASK | FLAG_GEODESIC_MASK; + if (flags & (~validFlags)) { return E_OPTION_INVALID; } + return E_SUCCESS; } diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index a5bf4c8510..8ca93906ce 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -235,7 +235,15 @@ It produces a collection of cells that are contained within the polygon. This function differs from `polygonToCells` in that it uses an experimental new algorithm which supports center-based, fully-contained, and -overlapping containment modes. +overlapping containment modes. You can additionally request a geodesic +traversal (great-circle edges) for very large polygons by enabling the +geodesic flag; this mode is significantly slower than the planar +algorithm, so prefer lower resolutions when using it for global-scale +objects. + +When geodesic traversal is enabled, only `CENTER`, `FULL`, and +`OVERLAPPING` containment modes are valid. `OVERLAPPING_BBOX` is invalid in +geodesic mode. @@ -289,6 +304,11 @@ The valid values for `flags` are: | `PolygonToCellsFlags.containment_overlapping` | Cell overlaps the shape at any point | `PolygonToCellsFlags.containment_overlapping_bbox` | Cell bounding box overlaps shape +When geodesic mode is enabled in Java, only +`PolygonToCellsFlags.containment_center`, +`PolygonToCellsFlags.containment_full`, and +`PolygonToCellsFlags.containment_overlapping` are valid. + @@ -317,6 +337,11 @@ The valid values for `flags` are: | `POLYGON_TO_CELLS_FLAGS.containment_overlapping` | Cell overlaps the shape at any point | `POLYGON_TO_CELLS_FLAGS.containment_overlapping_bbox` | Cell bounding box overlaps shape +When geodesic mode is enabled in JavaScript, only +`POLYGON_TO_CELLS_FLAGS.containment_center`, +`POLYGON_TO_CELLS_FLAGS.containment_full`, and +`POLYGON_TO_CELLS_FLAGS.containment_overlapping` are valid. + @@ -337,6 +362,8 @@ The valid values for `contain` are: | `overlap` | Cell overlaps the shape at any point | `bbox_overlap` | Cell bounding box overlaps shape +When geodesic mode is enabled in Python, only `center`, `full`, and `overlap` are valid. + @@ -353,6 +380,9 @@ The valid values for `containmentMode` are: | `h3.ContainmentOverlapping` | Cell overlaps the shape at any point | `h3.ContainmentOverlappingBbox` | Cell bounding box overlaps shape +When geodesic mode is enabled in Go, only `h3.ContainmentCenter`, +`h3.ContainmentFull`, and `h3.ContainmentOverlapping` are valid. + @@ -378,6 +408,8 @@ The valid values for the third `contain` argument are: | `overlap` | Cell overlaps the shape at any point | `bbox_overlap` | Cell bounding box overlaps shape +When geodesic mode is enabled in DuckDB, only `center`, `full`, and `overlap` are valid. + @@ -661,5 +693,3 @@ This function exists for memory management and is not exposed. - -