Geodesic polygon coverage #1052

holoskii · 2025-10-05T09:04:39Z

Summary

Adds the geodesic traversal mode to polygonToCellsExperimental, letting us follow true great-circle edges when filling massive polygons (think continental footprints and polar caps). The PR builds the geodesic stack from the ground up—new vector math helpers, bounding caps, iterator wiring, and plenty of validation coverage—while keeping the default planar behavior unchanged (and safely rejecting geodesic flags on the legacy API).

Highlights

New geodesic mode
- Introduces FLAG_GEODESIC_MASK (with helpers) so callers can opt into geodesic coverage; validates acceptable containment modes (FULL, OVERLAPPING) at both maxPolygonToCellsSizeExperimental and polygonToCellsExperimental.
- Dispatches a dedicated geodesic iterator path that reuses the existing traversal logic but evaluates cells using great-circle geometry.
Core geometry plumbing
- Adds GeodesicPolygon acceleration structures + helpers (geodesicPolygonCreate/Destroy, containment, boundary/cap intersection tests).
- Extends the vector toolkit (latLngToVec3, cross/dot/normalize, distance helpers) and new conversions: vec3dToCell, cellToVec3, cellToGeodesicBoundary, cellToSphereCap.
- Pulls bounding-cap tables into sphereCapTables.h for both prod code and tests; exposes SphereCap/AABB helpers.
Iterator & polyfill updates
- Refactors polygon iterators into polyfill_iterator.h (with _extra context) and wires in geodesic_iterator.c.
- Ensures geodesic flags short-circuit the legacy planar polygonToCells path.
Benchmarks, fuzzers, and tests
- Benchmarks: geodesic variants for existing benchmark.
- Fuzzers now exercise geodesic flag permutations.
- Added new test suites: testGeodesicPolygonToCellsExperimental, testSphereCap, testVec3d, plus extended testBBoxInternal, testPolygonToCells(Experimental) flag validation, and memory-failure coverage.
- Hooks new tests into CMakeLists.txt/CMakeTests.cmake.
Documentation
- Updates the Regions API docs to explain the geodesic flag, trade-offs, and usage tips.

Attribution

This geodesic polyfill work originated at Yellowbrick Data (where I’m currently employed). We’re upstreaming it to the H3 project so the broader community can benefit and use it for coverage.

CLAassistant · 2025-10-05T09:04:46Z

All committers have signed the CLA.

coveralls · 2025-10-06T17:52:37Z

coverage: 97.09% (-1.9%) from 98.946%
when pulling f3df47e on holoskii:geodesic_coverage
into 0ac92fb on uber:master.

ajfriend · 2025-10-09T01:02:26Z

Thanks @holoskii for the contribution! Geodesic polyfill is an exciting addition that I think a lot of folks would use, myself included!

Please bear with us as we try to find the time to properly review this PR. Some initial questions from my end:

Is it fair to characterize this PR as providing two major improvements:
1. geodesic mode for "polyfill", for polygons that would otherwise already be well-handled in the existing planar mode in polygonToCellsExperimental, and
2. the ability to do polyfill on polygons larger or tricker (e.g. around poles or the antimeridian) than those that could be handled in the previous algorithm (and, if so, would these improvements also generalize to the existing planar mode?)
This may become obvious after I'm able to review more closely, but how much code can we share between the planar and geodesic algorithms?
Is it necessary to expose Vec3d and GeodesicCellBoundary in the public api (h3api.h.in), or can we keep those structs internal?
You mentioned that this mode is slower than the existing algorithm. How much slower? Do you have any example benchmarks you could share?
Any thoughts on what to do about the failing CI checks? Or would you like some input from the team?

isaacbrodsky

Thanks for your contribution and opening this PR! Please let me know if you have questions. I am beginning to review the contribution and will add more comments as I read & test.

src/h3lib/include/h3api.h.in

isaacbrodsky · 2025-10-10T04:50:25Z

src/h3lib/lib/faceijk.c

+void _geoToHex2d(const LatLng *g, int res, int *face, Vec2d *v) {
+    Vec3d p = latLngToVec3(g);
+    _vec3dToHex2d(&p, res, face, v);
+}


From a code change perspective, I note this is a little surprising because it appears it affects the indexing code path as well. (Not necessarily wrong, just something to review in detail)

Extracted common bits into a helper to remove duplication. The original planar path behavior should be unchanged

isaacbrodsky · 2025-10-10T04:56:52Z

src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c

+        if (flags == CONTAINMENT_FULL || flags == CONTAINMENT_OVERLAPPING) {
+            uint32_t geodesicFlags = flags;
+            FLAG_SET_GEODESIC(geodesicFlags);
+            run(&geoPolygon, geodesicFlags, res);
+        }


We may want geodesic mode to be in a separate fuzzer to avoid fuzzer timeout issues.

Come to think of it, if we have more complicated input for the flags, we may want to also want to begin having a fuzzer for polygons (with or without holes) that also tries various flag values.

The fuzzer timeout was due to a bug, not algorithmic slowness. After the fix, there should be no more timeouts.

I kept planar and geodesic in a single fuzzer since they share the same entry point/interface. If runtime becomes an issue, we can first raise the timeout before splitting targets.

For flags, there’s a test that enumerates all combinations: TEST(invalidFlagsExperimental).

src/apps/testapps/testGeodesicPolygonInternal.c

src/h3lib/include/sphereCapTables.h

src/h3lib/include/polyfill_iterator.h

holoskii · 2025-10-24T09:03:57Z

Fixes / changes

Introduce an explicit epsilon for the SphereCap check. The previous iteration relied on tolerance built into precomputed cosine values; in practice, we also need an additional epsilon.
Use GeodesicPolygon* instead of void* in the iterator structure.
Fix the MSVC build by using a macro instead of a constant.
Move GeodesicCellBoundary and Vec3d to their own headers instead of defining them in the public API header.
Move internal headers to the include directory.
Add a new geodesic benchmark.
Rebase on the latest master.

Performance

I added a dedicated geodesic benchmark that exercises small, medium, and large shapes.

Small cells: up to 60× slower.
Medium cells / large shapes: ~20×.
“Primary” geodesic use case (huge shapes + large cells, e.g., NY ↔ London great-circle path): ~2.5×.

Notes:

At Yellowbrick Data, we’ve seen ~10× performance improvements by pushing heavy inlining (by moving code to headers), branch hints, etc. I avoided invasive changes here to keep code style intact; the current PR aims for minimal intrusion.

holoskii · 2025-10-24T09:06:43Z

Thanks @holoskii for the contribution! Geodesic polyfill is an exciting addition that I think a lot of folks would use, myself included!

Please bear with us as we try to find the time to properly review this PR. Some initial questions from my end:

* Is it fair to characterize this PR as providing two major improvements:
  
  1. geodesic mode for "polyfill", for polygons that would otherwise already be well-handled in the existing planar mode in `polygonToCellsExperimental`, and
  2. the ability to do polyfill on polygons larger or tricker (e.g. around poles or the antimeridian) than those that could be handled in the previous algorithm (and, if so, would these improvements also generalize to the existing planar mode?)

* This may become obvious after I'm able to review more closely, but how much code can we share between the planar and geodesic algorithms?

* Is it necessary to expose `Vec3d` and `GeodesicCellBoundary` in the public api (`h3api.h.in`), or can we keep those structs internal?

* You mentioned that this mode is slower than the existing algorithm. How much slower? Do you have any example benchmarks you could share?

* Any thoughts on what to do about the failing CI checks? Or would you like some input from the team?

1.1 - Yes.

1.2 - The geodesic approach behaves better for large polygons and near poles/antimeridian due to great-circle behavior, but it’s architecturally different from the planar path; there isn’t a practical feature we can “back-port” without re-implementing substantial pieces.

2 - We only share FaceIjk and iterator code between the two implementations. Intersection/containment logic and related optimization structures are written from scratch for the geodesic mode.

3 - Performance metrics are now included in the main PR comment (see “Performance” above) with the new benchmark.

holoskii · 2025-10-31T09:04:31Z

A few small fixes for the CI:

examples did not compile because of "constants.h" included in h3api.h. Removed that include to keep api header standalone
fuzzer hit a legitimate timeout. Fixed this by limiting resolution and number of points for geodesic runs

holoskii · 2025-11-03T21:45:07Z

@isaacbrodsky was right, separate fuzzer is a way to go. Old fuzzers left not enough time until timeout in each iteration. In the last commit I've added a separate geodesic fuzzer. PR is now fully ready for review

holoskii · 2025-11-12T12:04:47Z

Hi @nrabinowitz @dfellis @isaacbrodsky - friendly bump on this.
Current state: fixed CI issues, added benchmarks, answered questions and addressed the comments.

isaacbrodsky · 2025-11-14T01:01:28Z

Hi @nrabinowitz @dfellis @isaacbrodsky - friendly bump on this. Current state: fixed CI issues, added benchmarks, answered questions and addressed the comments.

Thanks for updating the PR and for the bump! We appreciate the contribution and know we need to review this, so it's just a matter of finding time to read through it in depth. We would like to bring this in for the next release of H3.

isaacbrodsky · 2025-12-11T20:21:56Z

CMakeLists.txt

+    src/h3lib/include/polyfill_iterator.h
+    src/h3lib/include/geodesic_iterator.h
+    src/h3lib/include/geodesic_polygon.h
+    src/h3lib/include/geodesic_polygon_internal.h
+    src/h3lib/include/geodesic_cell_boundary.h


very minor, but we can rename these files to camelCase as with other source files

isaacbrodsky · 2025-12-11T20:25:59Z

src/apps/benchmarks/benchmarkPolygonToCellsExperimentalGeodesic.c

+LatLng stateVerts[] = {{0.645778804554910, -1.903190792178713},
+                       {0.645682811446050, -1.780975856637062},
+                       {0.715595465293187, -1.781167842854781},
+                       {0.715578012000667, -1.903262350678044},
+                       {0.645778804554910, -1.903190792178713}};


I think either these coordinates are inverted compared to the London/NYC coordinates below, or other way around?

isaacbrodsky · 2025-12-11T20:29:33Z

src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c

+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;
+}


We could make a follow up issue to deduplicate the populateGeoLoop function between the three files its now in

isaacbrodsky · 2025-12-11T20:31:38Z

src/apps/fuzzers/fuzzerPolygonToCellsExperimentalGeodesic.c

+    if (res == 0) {
+        res = 1;  // resolution 1 tests more code paths compared to 0
+    }


If res=0 is a plausible input we should exercise it under the fuzzer too.

isaacbrodsky · 2025-12-11T20:43:39Z

src/apps/testapps/testH3Memory.c

+        uint32_t flags = CONTAINMENT_OVERLAPPING;
+        FLAG_SET_GEODESIC(flags);
+
+        resetMemoryCounters(1);


Does this need to test the case that maxPolygonToCellsSizeExperimental fails due to memory allocation failing?

isaacbrodsky · 2025-12-11T21:32:39Z

src/h3lib/lib/faceijk.c

+    // of a distortion vertex on the last edge
+    int additionalIteration = length == NUM_HEX_VERTS ? 1 : 0;
+
+    // convert each vertex to lat/lng


Suggested change

// convert each vertex to lat/lng

// convert each vertex to xyz

isaacbrodsky · 2025-12-11T21:35:14Z

src/h3lib/lib/faceijk.c

+ * Encodes a Vec3d coordinate to the corresponding icosahedral face and
+ * squared euclidean distance to that face center.
 *
- * @param g The spherical coordinates to encode.
- * @param res The desired H3 resolution for the encoding.
- * @param h The FaceIJK address of the containing cell at resolution res.
+ * @param v3d The Vec3d coordinates to encode.
+ * @param face The icosahedral face containing the coordinates.
+ * @param sqd The squared euclidean distance to its icosahedral face center.
 */
-void _geoToFaceIjk(const LatLng *g, int res, FaceIJK *h) {
-    // first convert to hex2d
-    Vec2d v;
-    _geoToHex2d(g, res, &h->face, &v);
+static void _vec3dToClosestFace(const Vec3d *v3d, int *face, double *sqd) {
+    // determine the icosahedron face
+    *face = 0;
+    // The distance between two farthest points is 2.0, therefore the square of
+    // the distance between two points should always be less or equal than 4.0 .
+    *sqd = 5.0;
+    for (int f = 0; f < NUM_ICOSA_FACES; ++f) {
+        double sqdT = vec3DistSq(&faceCenterPoint[f], v3d);
+        if (sqdT < *sqd) {
+            *face = f;
+            *sqd = sqdT;
+        }
+    }
+}


I would move this function to the bottom of the file if possible to minimize the diff

isaacbrodsky · 2025-12-11T21:42:53Z

src/h3lib/lib/geodesic_iterator.c

+        if (iter->_polygon->geoloop.numVerts <= 0) {
+            iterErrorPolygonCompact(iter, E_DOMAIN);
+            return NULL;
+        }


a zero-vertex geodesic polygon is not allowed?

isaacbrodsky · 2025-12-11T21:44:55Z

src/h3lib/lib/geodesic_iterator.c

+        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 <= 0) {
+                iterErrorPolygonCompact(iter, E_DOMAIN);
+                return NULL;
+            }
+        }
+        GeodesicPolygon *poly = geodesicPolygonCreate(iter->_polygon);
+        if (!poly) {
+            iterErrorPolygonCompact(iter, E_MEMORY_ALLOC);
+            return NULL;
+        }
+        iter->geodesicPoly = poly;


A lot of these cases don't seem to be tested. numHoles being set while holes is NULL makes sense to check although if someone mismatched the holse and numHoles there is not a lot we can do about it in general.

Regarding test coverage, do we want to update this PR to cover all cases, or merge this PR and then open additional PRs to complete the coverage of new code?

isaacbrodsky · 2025-12-11T21:52:30Z

src/h3lib/lib/geodesic_polygon.c

+        return E_DOMAIN;
+    }
+
+    *out = (GeodesicLoop){0};


Suggested change

*out = (GeodesicLoop){0};

*out = {0};

I assume this is not needed?

ajfriend requested review from dfellis, isaacbrodsky and nrabinowitz October 8, 2025 23:15

isaacbrodsky reviewed Oct 10, 2025

View reviewed changes

ajfriend added this to the v4.5 milestone Oct 16, 2025

isaacbrodsky reviewed Oct 16, 2025

View reviewed changes

holoskii added 11 commits October 24, 2025 10:03

Geodesic coverage

3759994

AUDIT_SOURCE_FILE_LIST

c041c98

Remove H3_EXPORT where not needed

881aecb

Move internal header to include dir

464bdbd

GeodesicCellBoundary to separate header

0e6c108

Attempt to fix MSVC compilation

fab4ab6

GeodesicPolygon * instead of void *

6d2acdf

Introduce epsilon for SphereCap check

db50026

sphereCapTables.h to include dir

0afa5b1

Add new benchmark

1fd3074

Audit src file list

11cf632

holoskii force-pushed the geodesic_coverage branch from c9be803 to 11cf632 Compare October 24, 2025 08:58

holoskii requested a review from isaacbrodsky October 27, 2025 09:48

holoskii added 3 commits October 29, 2025 10:50

Fix compilation: move helper function outside benchmark

685849e

Keep h3api.h free of non-standart includes

d3f0cc0

Slim down geodesic fuzzer

c7a24cd

Separate geodesic fuzzer

f3df47e

isaacbrodsky mentioned this pull request Nov 14, 2025

OGC-DGGRS-API feedback. New function getGeometricSubZones #1095

Open

isaacbrodsky reviewed Dec 11, 2025

View reviewed changes

	// convert each vertex to lat/lng
	// convert each vertex to xyz

Geodesic polygon coverage #1052

Are you sure you want to change the base?

Geodesic polygon coverage #1052

Conversation

holoskii commented Oct 5, 2025

Summary

Highlights

Attribution

Uh oh!

CLAassistant commented Oct 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

coveralls commented Oct 6, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ajfriend commented Oct 9, 2025

Uh oh!

isaacbrodsky left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

holoskii commented Oct 24, 2025

Fixes / changes

Performance

Uh oh!

holoskii commented Oct 24, 2025

Uh oh!

holoskii commented Oct 31, 2025

Uh oh!

holoskii commented Nov 3, 2025

Uh oh!

holoskii commented Nov 12, 2025

Uh oh!

isaacbrodsky commented Nov 14, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

CLAassistant commented Oct 5, 2025 •

edited

Loading

coveralls commented Oct 6, 2025 •

edited

Loading

isaacbrodsky left a comment •

edited

Loading