Vec3d indexing refactor

[ajfriend]

Summary

Vec3d refactor

Based on @holoskii 's work in #1052, we cherry pick out just the Vec3d indexing refactor.

The indexing path on master currently operates on "geo" lat/lng coordinates:

latLngToCell (old, on master)
├── _geoToFaceIjk → FaceIJK
│   ├── _geoToHex2d → face + Vec2d
│   │   ├── _geoToClosestFace (uses Vec3d internally, but discards it)
│   │   ├── _geoAzimuthRads
│   │   ├── gnomonic projection
│   │   └── polar → cartesian
│   └── _hex2dToCoordIJK → CoordIJK
└── _faceIjkToH3 → H3Index

The new path in this PR converts from lat/lng to Vec3d and does more operations in that representation:

vec3ToCell (new)
├── _vec3ToFaceIjk → FaceIJK
│   ├── _vec3ToHex2d → face + Vec2d
│   │   ├── _vec3ToClosestFace
│   │   ├── _vec3AzimuthRads
│   │   ├── gnomonic projection
│   │   └── polar → cartesian
│   └── _hex2dToCoordIJK → CoordIJK
└── _faceIjkToH3 → H3Index

latLngToCell (new, thin wrapper)
├── latLngToVec3
└── vec3ToCell

Additional changes

In addition to what was done in #1052, this PR makes everything go through the new Vec3d logic, for a complete migration:

• latLngToCell
• cellToLatLng
• cellToBoundary
• directedEdgeToBoundary
• vertexToLatLng

Moving everything let us delete all the "to geo" and "geo to" code, e.g., _faceIjkToGeo.

Also:

• Making a header-only vec3d.h provided a nice win: latLngToCell went from ~11% slower to ~3% faster than master
• 100% line and branch coverage on change

Cell assignment tests

Compares this PR to what's on master.

Code in src/apps/benchmarks/dumpCoreApi.c and stress.py (will be removed before landing). Run with uv run stress.py

Note: These are on top of the existing tests, bc*centers.txt and rand*centers.txt in tests/inputfiles/.

latLngToCell: random points

10,000 random points (uniform on sphere, fixed seed) tested at 8 resolutions (0, 1, 2, 3, 4, 7, 10, 15) = 80,000 tests.

Result: 0 mismatches

latLngToCell: boundary points

All cell boundary vertices and edge midpoints at resolutions 0–3, tested at 8 resolutions = 4,654,080 tests.

Result: ~12% mismatches. But this is expected, since we arbitrarily break ties between cells, and slightly different floating point math may break ties differently.

cellToLatLng and cellToBoundary

All cells at resolutions 0--4.

Result: all agree to within ~0.001 mm

Benchmarks

Code in src/apps/benchmarks/benchmarkCoreApi.c and bench.py. Run with uv run bench.py. (Will be removed before landing.)

Function	master	vec3-core	Change
latLngToCell	0.198 µs	0.192 µs	-3.4%
cellToLatLng	0.132 µs	0.109 µs	-17.6%
cellToBoundary	0.590 µs	0.443 µs	-25.0%
directedEdgeToBoundary	0.303 µs	0.258 µs	-14.7%
vertexToLatLng	0.172 µs	0.149 µs	-13.7%

Note that one of the bottlenecks of the new cellsToMultiPolygon code is directedEdgeToBoundary, so speedups here will help that function as well!

Follow-ups

Renaming

To help me understand the code, I had renamed a bunch of things, including using Vec2 and Vec3 consistently ("hex2d" and "vec2d" were used interchangeably, and function names didn't always match struct names). While that helped in working on the PR, it made this diff hard to read. I'll revisit in the future, since I think better naming would be nice if we continue to work on this code.

Header-only allows for helper function optimizations

Compiler optimizations across translation units makes a difference. By moving vec3d.h to be header-only and static inline, the latLngToCell benchmarks go from +11% (slower) to -3.4% (faster), compared to master. We can do the same elsewhere: coordijk.c has lots of small functions that get called in other .c files (translation units). I tested moving coordijk.c to be header-only, and got latLngToCell to be -9.7% (faster) compared to master, and cellToLatLng and cellToBoundary were both -33% (faster). We should do this in a follow-up PR. And there might be other files worth considering, like vec2d.h.

Trig-free rotation math

The rotation math in _vec3ToVec2 translates from 3d to angle and back to 3d. That rotation can be done without any trig functions. Instead of storing angles in faceAxesAzRadsCII (and note that we store all three, but only use the first!), we could instead store the Vec3 i-axis for each icosahedron face. This removes 5 trig function calls per indexing operation. Bigger change, separate PR.

Note, we can also remove some trig from this path in _vec3ToHex2d:

double r = acos(1 - sqd * 0.5);
r = tan(r);

by writing

double r = sqrt(sqd * (1.0 - sqd * 0.25)) / (1.0 - sqd * 0.5);

I tried this, but didn't see much movement in the benchmarks.

[coveralls]

coverage: 99.161% (+0.1%) from 99.066%
when pulling a06d168 on ajfriend:vec3d-core
into 006ba2c on uber:master.

[ajfriend]

Slightly different floating point math on Mac vs Linux. Reducing to 8 digits gets tests to pass on both platforms. (We do an exact string comparison in edgeLengthM.txt)