All committers have signed the CLA.

coverage: 98.543% (-0.6%) from 99.095%
when pulling fb857a1 on holoskii:geodesic_coverage
into 1d5346b on uber:master.

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?

Fixes / changes

1. 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.
2. Use GeodesicPolygon* instead of void* in the iterator structure.
3. Fix the MSVC build by using a macro instead of a constant.
4. Move GeodesicCellBoundary and Vec3d to their own headers instead of defining them in the public API header.
5. Move internal headers to the include directory.
6. Add a new geodesic benchmark.
7. 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.

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.

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

@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

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

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.

@holoskii Thanks for the update! I think CI is failing on the formatting assertion (needs make format with clang-format-14)

@holoskii Thanks for the update! I think CI is failing on the formatting assertion (needs make format with clang-format-14)

My update is a bit overdue, this year has been very busy so far 😅
I will address a few more comments (missed them in my first push), and then will use the correct formatter (I probably used wrong version)

Resolved final issues, formatted with version 14. Should be ready for the CI

2 small issues to raise after this PR is merged:

• deduplicate populateGeoLoop: Add geodesic polygon coverage mode for polygonToCellsExperimental #1052 (comment)
• add more modes to geodesic coverage Add geodesic polygon coverage mode for polygonToCellsExperimental #1052 (comment)

Windows workflow failed because on 32-bit Windows systems, the math library doesn't handle NaN inputs correctly. Instead of ignoring NaN data when building a bounding box, it lets the NaN poison the BBOX. This leads to us hitting NaN down the road and failing later.

I added a check at the very beginning, during polygon creation to catch invalid coordinates (NaN or Infinity). Now, we reject bad input right away with a clear error code (E_DOMAIN) instead of letting it break things further down the line.

Updated PR title and description, should be ready for the merge

Thanks again for this. This is awesome, and the effort is definitely appreciated.

I have a few questions about the interface, which don't necessarily need to block this PR, but I think we should discuss and resolve before releasing (or exiting "experimental" status). And since we definitely want to get this functionality over the line and released, we can also split this work up, so it isn't entirely on you, since you've already contributed quite a bit. :)

Discussion

I'm primarily interested in how this new mode will work with the "global" polygons (larger than a hemisphere, crossing poles and the antimeridian, etc.) that we can now output with cellsToMultiPolygon, and how "round trips" between these functions will behave.

Boundaries and center containment

For the purposes of testing (and because it'd be generally useful), would it be possible to add a CONTAINMENT_CENTER mode? (I'm hoping that this wouldn't be much work, and we could reuse much of what you already have here.) This would make the "round trip" testing easier, since the current modes require you to deal with boundary behavior. (Edit: I see it is mentioned here: #1052 (comment))

For example, if we take a single cell 0x8087fffffffffff, get its boundary, and run it through geodesic mode with
CONTAINMENT_FULL and CONTAINMENT_OVERLAPPING, we get 0 cells and 5 cells (the center cell and 4 of 6 neighboring cells), respectively. If we reverse the boundary, we get 0 cells and 6 cells (center and 5 of 6 neighbors).

Ideally, we'd get just the single cell back, but I understand boundary behavior makes that difficult with FULL and OVERLAPPING. (Side question: Given the known issues with boundaries, are these examples still behaving as you would expect?)

Polygon validation

Minimum number of points in a loop

Like I mention elsewhere in the PR, I think we should only accept polygons with loops containing at least 3 points, and return errors otherwise. There's no notion of "inside" or "outside" when you have only two points in a "polygon".

Winding order

For geodesic mode, I think we should require that polygons follow the "right-hand rule", with outer loops expected (and interpreted) to be in counter-clockwise winding order, and holes in clockwise order.

I recognize that this differs from what we do in the planar case, where we accept either ordering and figure out how to interpret it. But if we want to preserve the option of extending this mode to handle "global" polygons in the future, then being strict about loop winding is the only option. This would also align with the GeoJSON spec: RFC 7946, Section 3.1.6.

For example, if we take the boundary of a cell, reverse the ordering of the points so that they're in clockwise order, then this loop would represent "everything except the cell", and we'd want this function to interpret it that way. Currently, geodesic mode interprets the boundary and its reversal the same: as the smaller of the two options.

API guarantees and error codes

The public API docs should describe the supported input domain and what happens outside it. We should also test against that description.

For example, what polygons are guaranteed to work correctly? Which ones aren't? And, in the case of unsupported polygons, do we produce an error code, fail silently, or produce a partial result?

To start, I'd propose something like: "We support any polygon which is strictly contained within a hemisphere, and raise an error otherwise."

I think that can be implemented fairly easily in the current code by adding a check that all points within a polygon are in some half space, and the loops are oriented correctly (which we can do, for example, with the spherical area calculation code).

The plan for future versions of this function would be to relax the constraint on valid inputs, and eventually support all global polygons, like the ones we can get from cellsToMultiPolygon.

Testing

"Exhaustive" testing over the whole globe is something we've done with other functions that would be helpful here to ensure this is behaving as expected. For example, we could add some tests that evaluate all cells and their boundaries at, say, resolutions 0--4. We could also reverse the boundaries to check that we get the expected output or error code. We might also want to check grid rings and grid disks around the cells, for both small regions and regions approaching or exceeding a hemisphere.

In particular, I don't think we currently have any tests for what happens when the function receives a polygon larger than a hemisphere, do we?

Again, the CONTAINMENT_CENTER mode would be helpful here.

Next steps

That covers my thoughts and concerns on the current PR (which, again, is awesome. thanks!). It involves a few changes, but ones which I expect wouldn't involve a ton of work, since they're more about the API boundary than the function internals. Still, we should definitely have a discussion to see if folks agree, and then how we want to triage and distribute this work. It probably doesn't need to happen all in this PR.

Aside from my API comments above, I may have found a proper bug in the current implementation.

The test is the following: For a given cell h, I get the grid disk of some radius k around the cell, and convert it to a polygon with cellsToMultiPolygon(). I then run polygonToCellsExperimental() with CONTAINMENT_OVERLAPPING | FLAG_GEODESIC at the same resolution as the input cell. Recalling the issues mentioned above, we should expect that the output cells should be a superset of the grid disk.

That is the case for most cells, however, at res 0 and k=1, I'm seeing three test failures:

=== h = 0x8027fffffffffff,  k=1 ===
Input (7 cells):
  0x8027fffffffffff
  0x8013fffffffffff
  0x8029fffffffffff
  0x8049fffffffffff
  0x8045fffffffffff
  0x802bfffffffffff
  0x800ffffffffffff
Output (16 cells):
  0x8003fffffffffff
  0x8007fffffffffff
  0x800dfffffffffff
  0x800ffffffffffff
  0x8013fffffffffff
  0x801bfffffffffff
  0x801dfffffffffff
  0x8027fffffffffff
  0x8029fffffffffff
  0x802bfffffffffff
  0x8049fffffffffff
  0x804dfffffffffff
  0x8051fffffffffff
  0x8067fffffffffff
  0x806dfffffffffff
  0x806ffffffffffff

MISSING: 0x8045fffffffffff

=== h = 0x804bfffffffffff,  k=1 ===
Input (7 cells):
  0x8069fffffffffff
  0x804ffffffffffff
  0x802ffffffffffff
  0x8031fffffffffff
  0x8041fffffffffff
  0x8073fffffffffff
  0x804bfffffffffff
Output (13 cells):
  0x8017fffffffffff
  0x8025fffffffffff
  0x802ffffffffffff
  0x8033fffffffffff
  0x803dfffffffffff
  0x8041fffffffffff
  0x804bfffffffffff
  0x804ffffffffffff
  0x8065fffffffffff
  0x8069fffffffffff
  0x8073fffffffffff
  0x808dfffffffffff
  0x8095fffffffffff

MISSING: 0x8031fffffffffff

=== h = 0x80d3fffffffffff,  k=1 ===
Input (7 cells):
  0x80d3fffffffffff
  0x80e9fffffffffff
  0x80cffffffffffff
  0x80b7fffffffffff
  0x80b1fffffffffff
  0x80c7fffffffffff
  0x80e3fffffffffff
Output (15 cells):
  0x808ffffffffffff
  0x8093fffffffffff
  0x80a1fffffffffff
  0x80b1fffffffffff
  0x80b3fffffffffff
  0x80b7fffffffffff
  0x80c7fffffffffff
  0x80d3fffffffffff
  0x80d5fffffffffff
  0x80dffffffffffff
  0x80e3fffffffffff
  0x80e9fffffffffff
  0x80ebfffffffffff
  0x80effffffffffff
  0x80f3fffffffffff

MISSING: 0x80cffffffffffff

I can find more examples going up to res 1 and k = 2.

Edit: Here's the test code I was using: https://gist.github.com/ajfriend/88e366c5b464a15355d511633dcf5b73

And here's an example where just a single cell boundary (k=0) fails to recover the center cell:

=== h = 0x8081fffffffffff,  k=0 ===
Input (1 cells):
  0x8081fffffffffff
Output (4 cells):
  0x805ffffffffffff
  0x807dfffffffffff
  0x808bfffffffffff
  0x80a9fffffffffff

MISSING: 0x8081fffffffffff

{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "Polygon",
    "coordinates": [[
      [-44.140412918186037, 7.725284543519059],
      [-53.682044428130048, -0.170177724128589],
      [-51.776426167874277, -12.387611570633618],
      [-39.830275464745220, -16.779214626193742],
      [-30.225126243811296, -8.701739974658009],
      [-32.480128112549473, 3.454917535376012],
      [-44.140412918186037, 7.725284543519059]
    ]]
  }
}

Aside from my API comments above, I may have found a proper bug in the current implementation.

This is a legitimate issue!
It was caused by different paths used to get vertex vs edge value - results were differing by tiny amount, but that was enough to reject intersections. Fixed by using the same code path (even thought this is slightly slower). I have added a test that uses your approach of testing (cell -> disk -> polygon -> back to cells -> check).

@ajfriend #1052 (comment)
Great review! My reply to each of those:

Hemisphere approach

Now we actually try and find hemisphere that will fit out polygon and fail if we can't find one.
This hemisphere can actually be reused to fix another legitimate bug where for huge polygons the cell could "leak" to another hemisphere because of bad rejection approach. The new approach with hemisphere is more reliable.

CONTAINMENT_CENTER

Center containment mode was added to debug and fix issue with hemespheres. Also main iteration loop was optimized more, so we don't do expensive intersection test if we can avoid it

Boundaries issue

The behavior that you highlighted with different results for CW and CCW was not expected. This behavior occured because of the same issue causing wrong results: #1052 (comment). Fixed

Winding order

Even though I agree that having determenistic loop order would be a good change it would require doing a bigger refactor, and rethink our approach to current geodesic algorithm. As well as we would need to implement the same changes to the non-geodesic polyfill.

Additional testing

I have created several "full sweep" type of tests to try and find those nuanced bugs with FP errors.

• Add geodesic polygon coverage mode for polygonToCellsExperimental #1052 (comment)
• require full containment of cental cell within ring with k=1
• create hole in k=1 ring. Central cell with be excluded with full containment and included for partial

Thanks again for this work, addressing the comments, and fixing the issues, @holoskii. Really excited to get geodesic coverage into the library.

One note for alignment:

As well as we would need to implement the same changes to the non-geodesic polyfill

This isn't actually the case. We can't change the API of the non-geodesic polyfill (like adding a requirement for winding order) without a major version change. Since there are no plans for a major version bump any time soon, the non-geodesic mode is basically fixed as is. However, the geodesic mode is new, so we're free to design that API however we like. The point I was attempting to make is that we don't want paint ourselves into a similar corner with the API for geodesic mode. Since this code is tagged as experimental, we can be a bit more relaxed, but we'd want to avoid a final release that's agnostic to winding order if we're later planning on being strict about winding order to allow for global polygons.

Edit: And just to be clear, I don't think we need any changes for this PR, necessarily. I'm just calling this out for the future. For this PR, being agnostic to winding order is OK because 1) the function is still experimental, and 2) we could add a lightweight validation as a preprocessing step to check that polygons follow the right hand rule.

Just noting this here, since I was having trouble finding it earlier: #1052 (comment)

The updated PR now requires all loops to have >= 3 vertices. Thanks!

One chore I'd like to propose: splitting the faceijk.c / vec3d.c refactoring into its own PR, separate from the geodesic polygon code. It's a little more work to separate, but I think it's worth it for review quality and confidence. I'd be happy to help with this.

By splitting it up, we'll have an isolated PR focused on the internal refactor, which touches H3's core indexing functions. Subtle numerical changes here could potentially shift cell assignments, or affect indexing performance. @isaacbrodsky already flagged this in his review. Isolating it would let us add exhaustive round-trip tests and benchmarks to build full confidence, beyond what the existing test suite covers (since I'm not sure it was designed to capture regressions from this kind of change).

We can also write that PR with an eye towards reusable vec3d components for things like geodesic "line to cells" or "sphere cap to cells" functions.

It has the side benefit of making review easier. This is a significant (in size and functionality!) change, so splitting it into smaller PRs is helpful.

Curious to hear what other people think. And as mentioned, I'm happy to help out with this.

I'll review the recent changes, but my expectation at this point is that--pending this split--this PR is ready to merge.

The idea is great. This is a huge PR, so digesting it in smaller portions will indeed be easier. Thanks for the help with this part! I will be waiting for the focused faceijk/vec3 PR to review, and rebase this one after prep PR is merged

Great! I'll try to get tho this soon.

@holoskii, thanks for the patience: #1145

Please let me know if you have any suggestions.

PR has been dragging for too long. The new clean rebased PR will be opened