From acd2e4d9e79844362fd3cfc650ae8e0c2b2b42b6 Mon Sep 17 00:00:00 2001 From: Makar Ivashko <1212makar1212@gmail.com> Date: Mon, 22 Jun 2026 19:54:39 +0300 Subject: [PATCH 1/6] Geodesic polygon-to-cells coverage (squashed) Adds geodesic (great-circle) polygon-to-cells support to the polygonToCellsExperimental API via a FLAG_GEODESIC_MASK flag bit. Squashed from 39 commits on geodesic_coverage branch for clean rebase onto upstream/master. --- CHANGELOG.md | 1 + CMakeLists.txt | 21 +- CMakeTests.cmake | 7 + .../benchmarkPolygonToCellsExperimental.c | 3 + ...chmarkPolygonToCellsExperimentalGeodesic.c | 208 ++++++ src/apps/fuzzers/README.md | 2 +- .../fuzzerPolygonToCellsExperimental.c | 9 +- ...fuzzerPolygonToCellsExperimentalGeodesic.c | 121 ++++ .../fuzzerPolygonToCellsExperimentalNoHoles.c | 5 +- src/apps/testapps/testBBoxInternal.c | 56 ++ src/apps/testapps/testGeodesicCoverage.c | 109 ++++ .../testapps/testGeodesicGridDiskRoundtrip.c | 517 +++++++++++++++ .../testapps/testGeodesicPolygonInternal.c | 557 ++++++++++++++++ .../testGeodesicPolygonToCellsExperimental.c | 596 ++++++++++++++++++ src/apps/testapps/testGetIcosahedronFaces.c | 1 + src/apps/testapps/testH3IndexInternal.c | 46 ++ src/apps/testapps/testH3Memory.c | 261 +++++++- src/apps/testapps/testPolygonToCells.c | 4 + .../testapps/testPolygonToCellsExperimental.c | 100 ++- src/apps/testapps/testSphereCap.c | 199 ++++++ src/apps/testapps/testVec3d.c | 75 +++ src/h3lib/include/bbox.h | 44 ++ src/h3lib/include/faceijk.h | 5 + src/h3lib/include/geodesicCellBoundary.h | 35 + src/h3lib/include/geodesicIterator.h | 42 ++ src/h3lib/include/geodesicPolygonInternal.h | 108 ++++ src/h3lib/include/h3Index.h | 29 +- src/h3lib/include/h3api.h.in | 26 + src/h3lib/include/polyfill.h | 44 +- src/h3lib/include/polyfillIterator.h | 81 +++ src/h3lib/include/polygon.h | 4 - src/h3lib/include/sphereCapTables.h | 54 ++ src/h3lib/lib/algos.c | 8 + src/h3lib/lib/bbox.c | 79 +++ src/h3lib/lib/faceijk.c | 224 +++++++ src/h3lib/lib/geodesicIterator.c | 175 +++++ src/h3lib/lib/geodesicPolygon.c | 573 +++++++++++++++++ src/h3lib/lib/h3Index.c | 22 + src/h3lib/lib/polyfill.c | 47 +- src/h3lib/lib/polygon.c | 25 +- website/docs/api/regions.mdx | 34 +- 41 files changed, 4458 insertions(+), 99 deletions(-) create mode 100644 src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c create mode 100644 src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c create mode 100644 src/apps/testapps/testGeodesicCoverage.c create mode 100644 src/apps/testapps/testGeodesicGridDiskRoundtrip.c create mode 100644 src/apps/testapps/testGeodesicPolygonInternal.c create mode 100644 src/apps/testapps/testGeodesicPolygonToCellsExperimental.c create mode 100644 src/apps/testapps/testSphereCap.c create mode 100644 src/apps/testapps/testVec3d.c create mode 100644 src/h3lib/include/geodesicCellBoundary.h create mode 100644 src/h3lib/include/geodesicIterator.h create mode 100644 src/h3lib/include/geodesicPolygonInternal.h create mode 100644 src/h3lib/include/polyfillIterator.h create mode 100644 src/h3lib/include/sphereCapTables.h create mode 100644 src/h3lib/lib/geodesicIterator.c create mode 100644 src/h3lib/lib/geodesicPolygon.c 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..ea1a716286 --- /dev/null +++ b/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c @@ -0,0 +1,208 @@ +/* + * Copyright 2017, 2020-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @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..7b0b311fc0 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,11 @@ 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); + + // Don't run polygon without holes twice + 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..afba34ea5e --- /dev/null +++ b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c @@ -0,0 +1,121 @@ +/* + * Copyright 2023-2024 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..4b92fcdaa5 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) { @@ -441,4 +446,55 @@ 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"); + } } diff --git a/src/apps/testapps/testGeodesicCoverage.c b/src/apps/testapps/testGeodesicCoverage.c new file mode 100644 index 0000000000..ec1de4fe8e --- /dev/null +++ b/src/apps/testapps/testGeodesicCoverage.c @@ -0,0 +1,109 @@ +/* + * Copyright 2025 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..30fe5a3c57 --- /dev/null +++ b/src/apps/testapps/testGeodesicGridDiskRoundtrip.c @@ -0,0 +1,517 @@ +/* + * Copyright 2025 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..2821346d10 --- /dev/null +++ b/src/apps/testapps/testGeodesicPolygonInternal.c @@ -0,0 +1,557 @@ +/* + * Copyright 2025 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..ba195950fa --- /dev/null +++ b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c @@ -0,0 +1,596 @@ +/* + * Copyright 2025 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; +} + +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"); + } +} 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..7062bb668d --- /dev/null +++ b/src/apps/testapps/testVec3d.c @@ -0,0 +1,75 @@ +/* + * Copyright 2025 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..bc5eae0ef9 --- /dev/null +++ b/src/h3lib/include/geodesicCellBoundary.h @@ -0,0 +1,35 @@ +/* + * Copyright 2016-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 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..8b5b692aa8 --- /dev/null +++ b/src/h3lib/include/geodesicIterator.h @@ -0,0 +1,42 @@ +/* + * Copyright 2025 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..b65d9fe93f --- /dev/null +++ b/src/h3lib/include/geodesicPolygonInternal.h @@ -0,0 +1,108 @@ +/* + * Copyright 2025 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..8a3c68590c 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. @@ -92,47 +93,50 @@ /** * Gets the highest bit of the H3 index. */ -#define H3_GET_HIGH_BIT(h3) ((int)((((h3)&H3_HIGH_BIT_MASK) >> H3_MAX_OFFSET))) +#define H3_GET_HIGH_BIT(h3) \ + ((int)((((h3) & H3_HIGH_BIT_MASK) >> H3_MAX_OFFSET))) /** * Sets the highest bit of the h3 to v. */ -#define H3_SET_HIGH_BIT(h3, v) \ - (h3) = (((h3)&H3_HIGH_BIT_MASK_NEGATIVE) | \ +#define H3_SET_HIGH_BIT(h3, v) \ + (h3) = (((h3) & H3_HIGH_BIT_MASK_NEGATIVE) | \ (((uint64_t)(v)) << H3_MAX_OFFSET)) /** * Gets the integer mode of h3. */ -#define H3_GET_MODE(h3) ((int)((((h3)&H3_MODE_MASK) >> H3_MODE_OFFSET))) +#define H3_GET_MODE(h3) ((int)((((h3) & H3_MODE_MASK) >> H3_MODE_OFFSET))) /** * Sets the integer mode of h3 to v. */ #define H3_SET_MODE(h3, v) \ - (h3) = (((h3)&H3_MODE_MASK_NEGATIVE) | (((uint64_t)(v)) << H3_MODE_OFFSET)) + (h3) = \ + (((h3) & H3_MODE_MASK_NEGATIVE) | (((uint64_t)(v)) << H3_MODE_OFFSET)) /** * Gets the integer base cell of h3. */ -#define H3_GET_BASE_CELL(h3) ((int)((((h3)&H3_BC_MASK) >> H3_BC_OFFSET))) +#define H3_GET_BASE_CELL(h3) ((int)((((h3) & H3_BC_MASK) >> H3_BC_OFFSET))) /** * Sets the integer base cell of h3 to bc. */ #define H3_SET_BASE_CELL(h3, bc) \ - (h3) = (((h3)&H3_BC_MASK_NEGATIVE) | (((uint64_t)(bc)) << H3_BC_OFFSET)) + (h3) = (((h3) & H3_BC_MASK_NEGATIVE) | (((uint64_t)(bc)) << H3_BC_OFFSET)) /** * Gets the integer resolution of h3. */ -#define H3_GET_RESOLUTION(h3) ((int)((((h3)&H3_RES_MASK) >> H3_RES_OFFSET))) +#define H3_GET_RESOLUTION(h3) ((int)((((h3) & H3_RES_MASK) >> H3_RES_OFFSET))) /** * Sets the integer resolution of h3. */ #define H3_SET_RESOLUTION(h3, res) \ - (h3) = (((h3)&H3_RES_MASK_NEGATIVE) | (((uint64_t)(res)) << H3_RES_OFFSET)) + (h3) = \ + (((h3) & H3_RES_MASK_NEGATIVE) | (((uint64_t)(res)) << H3_RES_OFFSET)) /** * Gets the resolution res integer digit (0-7) of h3. @@ -145,15 +149,15 @@ * Sets a value in the reserved space. Setting to non-zero may produce invalid * indexes. */ -#define H3_SET_RESERVED_BITS(h3, v) \ - (h3) = (((h3)&H3_RESERVED_MASK_NEGATIVE) | \ +#define H3_SET_RESERVED_BITS(h3, v) \ + (h3) = (((h3) & H3_RESERVED_MASK_NEGATIVE) | \ (((uint64_t)(v)) << H3_RESERVED_OFFSET)) /** * Gets a value in the reserved space. Should always be zero for valid indexes. */ #define H3_GET_RESERVED_BITS(h3) \ - ((int)((((h3)&H3_RESERVED_MASK) >> H3_RESERVED_OFFSET))) + ((int)((((h3) & H3_RESERVED_MASK) >> H3_RESERVED_OFFSET))) /** * Sets the resolution res digit of h3 to the integer digit (0-7) @@ -181,5 +185,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..335c11f167 100644 --- a/src/h3lib/include/polyfill.h +++ b/src/h3lib/include/polyfill.h @@ -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..e6f2577669 --- /dev/null +++ b/src/h3lib/include/polyfillIterator.h @@ -0,0 +1,81 @@ +/* + * Copyright 2023 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file 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..73a65ad9e5 --- /dev/null +++ b/src/h3lib/include/sphereCapTables.h @@ -0,0 +1,54 @@ +/* + * 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 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..ff2f5dfe67 100644 --- a/src/h3lib/lib/bbox.c +++ b/src/h3lib/lib/bbox.c @@ -23,9 +23,12 @@ #include #include +#include "assert.h" #include "constants.h" #include "h3Index.h" #include "latLng.h" +#include "sphereCapTables.h" +#include "vec3d.h" /** * Width of the bounding box, in rads @@ -309,3 +312,79 @@ 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]; + + // Find the vector in the great circle plane that is perpendicular + // to the current Cartesian axis. This vector points to a potential + // extremum. + Vec3d extremum_dir = vec3Cross(*axis, normal); + 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..1697a17ed0 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; + 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..ffd3c8726d --- /dev/null +++ b/src/h3lib/lib/geodesicIterator.c @@ -0,0 +1,175 @@ +/* + * Copyright 2025 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) { + // Center is inside (guaranteed above). Fully inside only if + // boundaries do not intersect + if (!boundaryIntersection) { + iter->cell = cell; + return; + } + cell = nextCell(cell); + continue; + } + + // 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..b05ae81a45 --- /dev/null +++ b/src/h3lib/lib/geodesicPolygon.c @@ -0,0 +1,573 @@ +/* + * Copyright 2025 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" +#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. Use 1e-28 — 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. + const double sideEps = 1e-28; + if ((b1Side * b2Side > sideEps) || (a1Side * a2Side > sideEps)) { + 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..1ca00cc4a0 100644 --- a/src/h3lib/lib/polyfill.c +++ b/src/h3lib/lib/polyfill.c @@ -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 @@ -44,12 +46,6 @@ * 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 +273,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 +282,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 +327,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 +341,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 +427,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 +603,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 +753,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 +771,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..32a8649e2c 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 `FULL` and `OVERLAPPING` +containment modes are valid. `CENTER` and `OVERLAPPING_BBOX` are invalid in +geodesic mode. @@ -289,6 +304,10 @@ 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_full` and +`PolygonToCellsFlags.containment_overlapping` are valid. + @@ -317,6 +336,10 @@ 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_full` and +`POLYGON_TO_CELLS_FLAGS.containment_overlapping` are valid. + @@ -337,6 +360,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 `full` and `overlap` are valid. + @@ -353,6 +378,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.ContainmentFull` and +`h3.ContainmentOverlapping` are valid. + @@ -378,6 +406,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 `full` and `overlap` are valid. + @@ -661,5 +691,3 @@ This function exists for memory management and is not exposed. - - From b8dd07188d391827c9eed88aac8485c22d84f59a Mon Sep 17 00:00:00 2001 From: Makar Ivashko <1212makar1212@gmail.com> Date: Mon, 22 Jun 2026 20:19:39 +0300 Subject: [PATCH 2/6] Housekeeping: clang-format, dead includes, stale comment, copyright years --- src/apps/testapps/testGeodesicCoverage.c | 2 +- .../testapps/testGeodesicGridDiskRoundtrip.c | 2 +- .../testapps/testGeodesicPolygonInternal.c | 2 +- .../testGeodesicPolygonToCellsExperimental.c | 2 +- src/apps/testapps/testVec3d.c | 2 +- src/h3lib/include/geodesicCellBoundary.h | 2 +- src/h3lib/include/geodesicIterator.h | 2 +- src/h3lib/include/geodesicPolygonInternal.h | 2 +- src/h3lib/include/h3Index.h | 27 +++++++++---------- src/h3lib/include/polyfill.h | 2 +- src/h3lib/include/polyfillIterator.h | 2 +- src/h3lib/include/sphereCapTables.h | 2 +- src/h3lib/lib/bbox.c | 3 +-- src/h3lib/lib/geodesicIterator.c | 2 +- src/h3lib/lib/geodesicPolygon.c | 2 +- src/h3lib/lib/polyfill.c | 7 +---- 16 files changed, 27 insertions(+), 36 deletions(-) diff --git a/src/apps/testapps/testGeodesicCoverage.c b/src/apps/testapps/testGeodesicCoverage.c index ec1de4fe8e..3f7f6a5c59 100644 --- a/src/apps/testapps/testGeodesicCoverage.c +++ b/src/apps/testapps/testGeodesicCoverage.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/apps/testapps/testGeodesicGridDiskRoundtrip.c b/src/apps/testapps/testGeodesicGridDiskRoundtrip.c index 30fe5a3c57..3feaa07026 100644 --- a/src/apps/testapps/testGeodesicGridDiskRoundtrip.c +++ b/src/apps/testapps/testGeodesicGridDiskRoundtrip.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/apps/testapps/testGeodesicPolygonInternal.c b/src/apps/testapps/testGeodesicPolygonInternal.c index 2821346d10..e64ae2ecd0 100644 --- a/src/apps/testapps/testGeodesicPolygonInternal.c +++ b/src/apps/testapps/testGeodesicPolygonInternal.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c index ba195950fa..d166653a93 100644 --- a/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c +++ b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/apps/testapps/testVec3d.c b/src/apps/testapps/testVec3d.c index 7062bb668d..e2aef5ffde 100644 --- a/src/apps/testapps/testVec3d.c +++ b/src/apps/testapps/testVec3d.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/include/geodesicCellBoundary.h b/src/h3lib/include/geodesicCellBoundary.h index bc5eae0ef9..6bf6458536 100644 --- a/src/h3lib/include/geodesicCellBoundary.h +++ b/src/h3lib/include/geodesicCellBoundary.h @@ -1,5 +1,5 @@ /* - * Copyright 2016-2021 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/include/geodesicIterator.h b/src/h3lib/include/geodesicIterator.h index 8b5b692aa8..17c818ec2d 100644 --- a/src/h3lib/include/geodesicIterator.h +++ b/src/h3lib/include/geodesicIterator.h @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/include/geodesicPolygonInternal.h b/src/h3lib/include/geodesicPolygonInternal.h index b65d9fe93f..d29f2bfd76 100644 --- a/src/h3lib/include/geodesicPolygonInternal.h +++ b/src/h3lib/include/geodesicPolygonInternal.h @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/include/h3Index.h b/src/h3lib/include/h3Index.h index 8a3c68590c..74773d8507 100644 --- a/src/h3lib/include/h3Index.h +++ b/src/h3lib/include/h3Index.h @@ -93,50 +93,47 @@ /** * Gets the highest bit of the H3 index. */ -#define H3_GET_HIGH_BIT(h3) \ - ((int)((((h3) & H3_HIGH_BIT_MASK) >> H3_MAX_OFFSET))) +#define H3_GET_HIGH_BIT(h3) ((int)((((h3)&H3_HIGH_BIT_MASK) >> H3_MAX_OFFSET))) /** * Sets the highest bit of the h3 to v. */ -#define H3_SET_HIGH_BIT(h3, v) \ - (h3) = (((h3) & H3_HIGH_BIT_MASK_NEGATIVE) | \ +#define H3_SET_HIGH_BIT(h3, v) \ + (h3) = (((h3)&H3_HIGH_BIT_MASK_NEGATIVE) | \ (((uint64_t)(v)) << H3_MAX_OFFSET)) /** * Gets the integer mode of h3. */ -#define H3_GET_MODE(h3) ((int)((((h3) & H3_MODE_MASK) >> H3_MODE_OFFSET))) +#define H3_GET_MODE(h3) ((int)((((h3)&H3_MODE_MASK) >> H3_MODE_OFFSET))) /** * Sets the integer mode of h3 to v. */ #define H3_SET_MODE(h3, v) \ - (h3) = \ - (((h3) & H3_MODE_MASK_NEGATIVE) | (((uint64_t)(v)) << H3_MODE_OFFSET)) + (h3) = (((h3)&H3_MODE_MASK_NEGATIVE) | (((uint64_t)(v)) << H3_MODE_OFFSET)) /** * Gets the integer base cell of h3. */ -#define H3_GET_BASE_CELL(h3) ((int)((((h3) & H3_BC_MASK) >> H3_BC_OFFSET))) +#define H3_GET_BASE_CELL(h3) ((int)((((h3)&H3_BC_MASK) >> H3_BC_OFFSET))) /** * Sets the integer base cell of h3 to bc. */ #define H3_SET_BASE_CELL(h3, bc) \ - (h3) = (((h3) & H3_BC_MASK_NEGATIVE) | (((uint64_t)(bc)) << H3_BC_OFFSET)) + (h3) = (((h3)&H3_BC_MASK_NEGATIVE) | (((uint64_t)(bc)) << H3_BC_OFFSET)) /** * Gets the integer resolution of h3. */ -#define H3_GET_RESOLUTION(h3) ((int)((((h3) & H3_RES_MASK) >> H3_RES_OFFSET))) +#define H3_GET_RESOLUTION(h3) ((int)((((h3)&H3_RES_MASK) >> H3_RES_OFFSET))) /** * Sets the integer resolution of h3. */ #define H3_SET_RESOLUTION(h3, res) \ - (h3) = \ - (((h3) & H3_RES_MASK_NEGATIVE) | (((uint64_t)(res)) << H3_RES_OFFSET)) + (h3) = (((h3)&H3_RES_MASK_NEGATIVE) | (((uint64_t)(res)) << H3_RES_OFFSET)) /** * Gets the resolution res integer digit (0-7) of h3. @@ -149,15 +146,15 @@ * Sets a value in the reserved space. Setting to non-zero may produce invalid * indexes. */ -#define H3_SET_RESERVED_BITS(h3, v) \ - (h3) = (((h3) & H3_RESERVED_MASK_NEGATIVE) | \ +#define H3_SET_RESERVED_BITS(h3, v) \ + (h3) = (((h3)&H3_RESERVED_MASK_NEGATIVE) | \ (((uint64_t)(v)) << H3_RESERVED_OFFSET)) /** * Gets a value in the reserved space. Should always be zero for valid indexes. */ #define H3_GET_RESERVED_BITS(h3) \ - ((int)((((h3) & H3_RESERVED_MASK) >> H3_RESERVED_OFFSET))) + ((int)((((h3)&H3_RESERVED_MASK) >> H3_RESERVED_OFFSET))) /** * Sets the resolution res digit of h3 to the integer digit (0-7) diff --git a/src/h3lib/include/polyfill.h b/src/h3lib/include/polyfill.h index 335c11f167..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. diff --git a/src/h3lib/include/polyfillIterator.h b/src/h3lib/include/polyfillIterator.h index e6f2577669..333c68f455 100644 --- a/src/h3lib/include/polyfillIterator.h +++ b/src/h3lib/include/polyfillIterator.h @@ -1,5 +1,5 @@ /* - * Copyright 2023 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/include/sphereCapTables.h b/src/h3lib/include/sphereCapTables.h index 73a65ad9e5..48da1054f9 100644 --- a/src/h3lib/include/sphereCapTables.h +++ b/src/h3lib/include/sphereCapTables.h @@ -1,5 +1,5 @@ /* - * Copyright 2016-2018, 2020-2021 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/lib/bbox.c b/src/h3lib/lib/bbox.c index ff2f5dfe67..cebe2cea8b 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. @@ -23,7 +23,6 @@ #include #include -#include "assert.h" #include "constants.h" #include "h3Index.h" #include "latLng.h" diff --git a/src/h3lib/lib/geodesicIterator.c b/src/h3lib/lib/geodesicIterator.c index ffd3c8726d..2bea03d2c4 100644 --- a/src/h3lib/lib/geodesicIterator.c +++ b/src/h3lib/lib/geodesicIterator.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/lib/geodesicPolygon.c b/src/h3lib/lib/geodesicPolygon.c index b05ae81a45..1aa40b57e3 100644 --- a/src/h3lib/lib/geodesicPolygon.c +++ b/src/h3lib/lib/geodesicPolygon.c @@ -1,5 +1,5 @@ /* - * Copyright 2025 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/lib/polyfill.c b/src/h3lib/lib/polyfill.c index 1ca00cc4a0..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. @@ -41,11 +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. - */ /** All cells that contain the north pole, by res */ static H3Index NORTH_POLE_CELLS[MAX_H3_RES + 1] = { 0x8001fffffffffff, 0x81033ffffffffff, 0x820327fffffffff, 0x830326fffffffff, From 41540394f2bee34a15b2ebe7f69c9b11d17875bc Mon Sep 17 00:00:00 2001 From: Makar Ivashko <1212makar1212@gmail.com> Date: Tue, 23 Jun 2026 11:43:24 +0300 Subject: [PATCH 3/6] Fix maybe-uninitialized lastFijk in _faceIjkPentToCellBoundaryGeodesic gcc -O2 inlines _ijkToHex2d and sees the uninitialized read. Initialize to zero like the non-pentagon counterpart already does. --- src/h3lib/lib/faceijk.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/h3lib/lib/faceijk.c b/src/h3lib/lib/faceijk.c index 1697a17ed0..a87af1a63f 100644 --- a/src/h3lib/lib/faceijk.c +++ b/src/h3lib/lib/faceijk.c @@ -652,7 +652,7 @@ void _faceIjkPentToCellBoundaryGeodesic(const FaceIJK *h, int res, int start, // adjust the face of each vertex as appropriate and introduce // edge-crossing vertices as needed g->numVerts = 0; - FaceIJK lastFijk; + FaceIJK lastFijk = {0}; for (int vert = start; vert < start + length + additionalIteration; vert++) { int v = vert % NUM_PENT_VERTS; From 163a12c727b51f9ebac42afd97f5c0ece48e469d Mon Sep 17 00:00:00 2001 From: Makar Ivashko <1212makar1212@gmail.com> Date: Tue, 23 Jun 2026 13:41:12 +0300 Subject: [PATCH 4/6] Fix docs: geodesic mode permits CENTER containment The docs stated that only FULL and OVERLAPPING are valid when the geodesic flag is set. This contradicts the public API, which permits CENTER as well: h3api.h.in documents it, validatePolygonFlags accepts it (polygon.c), geodesicIteratorStep implements it (geodesicIterator.c), and testGeodesicPolygonToCellsExperimental.c exercises it. Only OVERLAPPING_BBOX is invalid in geodesic mode. Update the umbrella statement and all per-binding sections (Java, JavaScript, Python, Go, DuckDB) to list CENTER as valid. --- website/docs/api/regions.mdx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 32a8649e2c..8ca93906ce 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -241,8 +241,8 @@ 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 `FULL` and `OVERLAPPING` -containment modes are valid. `CENTER` and `OVERLAPPING_BBOX` are invalid in +When geodesic traversal is enabled, only `CENTER`, `FULL`, and +`OVERLAPPING` containment modes are valid. `OVERLAPPING_BBOX` is invalid in geodesic mode. @@ -337,7 +338,8 @@ The valid values for `flags` are: | `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_full` and +`POLYGON_TO_CELLS_FLAGS.containment_center`, +`POLYGON_TO_CELLS_FLAGS.containment_full`, and `POLYGON_TO_CELLS_FLAGS.containment_overlapping` are valid. @@ -360,7 +362,7 @@ 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 `full` and `overlap` are valid. +When geodesic mode is enabled in Python, only `center`, `full`, and `overlap` are valid. @@ -378,8 +380,8 @@ 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.ContainmentFull` and -`h3.ContainmentOverlapping` are valid. +When geodesic mode is enabled in Go, only `h3.ContainmentCenter`, +`h3.ContainmentFull`, and `h3.ContainmentOverlapping` are valid. @@ -406,7 +408,7 @@ 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 `full` and `overlap` are valid. +When geodesic mode is enabled in DuckDB, only `center`, `full`, and `overlap` are valid. From 8cc4926cdd074ff4249665f6a54ca07745b7389d Mon Sep 17 00:00:00 2001 From: Makar Ivashko <1212makar1212@gmail.com> Date: Tue, 23 Jun 2026 13:59:25 +0300 Subject: [PATCH 5/6] Fix geodesic correctness: arc extrema + boundary vertex logic --- src/apps/testapps/testBBoxInternal.c | 38 ++++++++++ .../testGeodesicPolygonToCellsExperimental.c | 76 +++++++++++++++++++ src/h3lib/lib/bbox.c | 9 ++- src/h3lib/lib/geodesicIterator.c | 37 +++++++-- 4 files changed, 149 insertions(+), 11 deletions(-) diff --git a/src/apps/testapps/testBBoxInternal.c b/src/apps/testapps/testBBoxInternal.c index 4b92fcdaa5..8072c860ba 100644 --- a/src/apps/testapps/testBBoxInternal.c +++ b/src/apps/testapps/testBBoxInternal.c @@ -51,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}}; @@ -497,4 +504,35 @@ SUITE(BBox) { 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/testGeodesicPolygonToCellsExperimental.c b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c index d166653a93..252bbfb63e 100644 --- a/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c +++ b/src/apps/testapps/testGeodesicPolygonToCellsExperimental.c @@ -240,6 +240,14 @@ static int64_t _geodesicFillSorted(const GeoPolygon *polygon, int res, 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; @@ -593,4 +601,72 @@ SUITE(geodesicPolygonToCellsExperimental) { 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/h3lib/lib/bbox.c b/src/h3lib/lib/bbox.c index cebe2cea8b..9bf74cc067 100644 --- a/src/h3lib/lib/bbox.c +++ b/src/h3lib/lib/bbox.c @@ -356,10 +356,11 @@ void aabbUpdateWithArcExtrema(AABB *aabb, const Vec3d *v1, const Vec3d *v2, for (int i = 0; i < 3; ++i) { const Vec3d *axis = &UNIT_AXES[i]; - // Find the vector in the great circle plane that is perpendicular - // to the current Cartesian axis. This vector points to a potential - // extremum. - Vec3d extremum_dir = vec3Cross(*axis, normal); + // 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 diff --git a/src/h3lib/lib/geodesicIterator.c b/src/h3lib/lib/geodesicIterator.c index 2bea03d2c4..6c3b51c2b2 100644 --- a/src/h3lib/lib/geodesicIterator.c +++ b/src/h3lib/lib/geodesicIterator.c @@ -131,14 +131,37 @@ void geodesicIteratorStep(IterCellsPolygonCompact *iter, H3Index cell) { geodesicPolygonBoundaryIntersects(poly, &boundary, &cap); if (mode == CONTAINMENT_FULL) { - // Center is inside (guaranteed above). Fully inside only if - // boundaries do not intersect - if (!boundaryIntersection) { - iter->cell = cell; - return; + // 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; } - 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 From b7d298225f36c4e8de8e3514e096d9b2c6be9ac3 Mon Sep 17 00:00:00 2001 From: Makar Ivashko <1212makar1212@gmail.com> Date: Mon, 6 Jul 2026 11:53:01 +0300 Subject: [PATCH 6/6] Address PR #1178 review feedback from dfellis - Fix copyright years on new files (benchmark + fuzzer geodesic) - Clarify fuzzer duplicate-run skip comment - Promote sideEps to file-scope SIDE_EPS constant in geodesicPolygon.c Co-Authored-By: Claude Opus 4.6 (1M context) --- ...enchmarkPolygonToCellsExperimentalGeodesic.c | 2 +- .../fuzzers/fuzzerPolygonToCellsExperimental.c | 3 ++- .../fuzzerPolygonToCellsExperimentalGeodesic.c | 2 +- src/h3lib/lib/geodesicPolygon.c | 17 ++++++++++------- 4 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c b/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c index ea1a716286..2bbf89fa2b 100644 --- a/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c +++ b/src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c @@ -1,5 +1,5 @@ /* - * Copyright 2017, 2020-2021 Uber Technologies, Inc. + * 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. diff --git a/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c b/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c index 7b0b311fc0..a1302940f2 100644 --- a/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c +++ b/src/apps/fuzzers/fuzzerPolygonToCellsExperimental.c @@ -96,7 +96,8 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { geoPolygon.numHoles = originalNumHoles; run(&geoPolygon, flags, res); - // Don't run polygon without holes twice + // If polygon has no holes - it has already been tested + // Don't test it again if (originalNumHoles == 0) { continue; } diff --git a/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c index afba34ea5e..388f01f8f7 100644 --- a/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c +++ b/src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c @@ -1,5 +1,5 @@ /* - * Copyright 2023-2024 Uber Technologies, Inc. + * 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. diff --git a/src/h3lib/lib/geodesicPolygon.c b/src/h3lib/lib/geodesicPolygon.c index 1aa40b57e3..45a94d6008 100644 --- a/src/h3lib/lib/geodesicPolygon.c +++ b/src/h3lib/lib/geodesicPolygon.c @@ -22,6 +22,10 @@ #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" @@ -49,13 +53,12 @@ static bool _geodesicEdgesCross(const Vec3d *a1, const Vec3d *a2, // ~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. Use 1e-28 — 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. - const double sideEps = 1e-28; - if ((b1Side * b2Side > sideEps) || (a1Side * a2Side > sideEps)) { + // 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; }