diff --git a/CMakeLists.txt b/CMakeLists.txt index 9426b868bc..3ba251c032 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -184,6 +184,7 @@ set(LIB_SOURCE_FILES src/h3lib/lib/localij.c src/h3lib/lib/latLng.c src/h3lib/lib/directedEdge.c + src/h3lib/lib/iterGosper.c src/h3lib/lib/mathExtensions.c src/h3lib/lib/iterators.c src/h3lib/lib/faceijk.c @@ -281,6 +282,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testGosperIter.c src/apps/testapps/testGeoLoopArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c @@ -320,6 +322,7 @@ set(OTHER_SOURCE_FILES src/apps/benchmarks/benchmarkGridDiskCells.c src/apps/benchmarks/benchmarkGridPathCells.c src/apps/benchmarks/benchmarkDirectedEdge.c + src/apps/benchmarks/benchmarkGosperIter.c src/apps/benchmarks/benchmarkVertex.c src/apps/benchmarks/benchmarkIsValidCell.c src/apps/benchmarks/benchmarkH3Api.c @@ -672,11 +675,17 @@ if(BUILD_BENCHMARKS) src/apps/benchmarks/benchmarkGridPathCells.c) add_h3_benchmark(benchmarkDirectedEdge src/apps/benchmarks/benchmarkDirectedEdge.c) + if(ENABLE_REQUIRES_ALL_SYMBOLS) + add_h3_benchmark(benchmarkGosperIter + src/apps/benchmarks/benchmarkGosperIter.c) + endif() add_h3_benchmark(benchmarkVertex src/apps/benchmarks/benchmarkVertex.c) add_h3_benchmark(benchmarkIsValidCell src/apps/benchmarks/benchmarkIsValidCell.c) - add_h3_benchmark(benchmarkCellsToPolyAlgos - src/apps/benchmarks/benchmarkCellsToPolyAlgos.c) + if(ENABLE_REQUIRES_ALL_SYMBOLS) + add_h3_benchmark(benchmarkCellsToPolyAlgos + src/apps/benchmarks/benchmarkCellsToPolyAlgos.c) + endif() add_h3_benchmark(benchmarkCellToChildren src/apps/benchmarks/benchmarkCellToChildren.c) add_h3_benchmark(benchmarkPolygonToCells diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 7500b47297..f48c8f9315 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -260,6 +260,7 @@ add_h3_test(testH3IteratorsInternal src/apps/testapps/testH3IteratorsInternal.c) add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) +add_h3_test(testGosperIter src/apps/testapps/testGosperIter.c) add_h3_test(testGeoLoopArea src/apps/testapps/testGeoLoopArea.c) add_h3_test_with_arg(testH3NeighborRotations diff --git a/branch_status.md b/branch_status.md new file mode 100644 index 0000000000..f6c8776c5e --- /dev/null +++ b/branch_status.md @@ -0,0 +1,42 @@ +# Branch: `child_edge_iter` + +**Date:** 2026-03-08 21:01 PDT + +## Summary + +This branch adds a Gosper island boundary iterator (`IterGosper`) and integrates it into the `cellsToMultiPolygon` pipeline. The iterator walks the fractal boundary of a cell's child set at a given resolution, yielding directed edges in geometric (CCW) order without enumerating any internal edges. + +## What was done + +### Gosper iterator + +- **`gosperIter.h` / `gosperIter.c`** — New iterator that yields the directed edges forming the outline of a cell's child set at a target resolution. Works for any parent/child resolution gap (0 to 15), hexagons and pentagons. Uses a boundary walk algorithm: at each resolution level, the outline of a parent's child set follows a repeating 18-step cycle (3 boundary segments per face, 6 faces). +- **`mathExtensions.h`** — Added `_ipow()` helper for integer exponentiation. +- **`testGosperIter.c`** — Exhaustive tests for res 0-4 (all cells x all child resolutions), plus high-res tests covering hex/pentagon at res 5/8/10 with gaps reaching res 15. Verifies edge count, validity, boundary membership, and loop connectivity. Includes expected-edge-sequence tests generated by an independent implementation. +- **`benchmarkGosperIter.c`** — Benchmarks for the iterator at various resolution gaps. +- **`CMakeLists.txt` / `CMakeTests.cmake`** — Build integration for tests and benchmarks. + +### Integration into cellsToMultiPolygon + +- **`cellsToMultiPoly.c`** — Added `cellsToMultiPolygonGosper()`, a new entry point that accepts possibly-compacted (mixed-resolution) cells and a target resolution. Uses the Gosper iterator to enumerate only boundary edges, then feeds them through the existing cancel-pairs / extract-loops / build-polygons pipeline. Extracted shared `arcSetToMultiPolygon()` helper to eliminate duplication between the flat and Gosper paths. Added overflow check for safety. +- **`cellsToMultiPoly.h`** — Declared `cellsToMultiPolygonGosper()`. +- **`benchmarkCellsToPolyAlgos.c`** — Added Colorado-based benchmarks at res 6, 7, and 8 comparing direct (flat), Gosper (flat), and Gosper (compacted) approaches. Extracted `coloradoCells()` and `compactAndCount()` helpers. + +## Benchmark results + +Compacted cells + Gosper iterator vs. flat `cellsToMultiPolygon` on Colorado: + +| Resolution | Direct (flat) | Gosper (flat) | Gosper (compact) | Compact speedup | +|------------|---------------|---------------|------------------|-----------------| +| 6 | 1,932 us | 1,711 us | 537 us | 3.6x | +| 7 | 24,371 us | 20,367 us | 1,908 us | 12.8x | +| 8 | 265,051 us | 240,309 us | 7,640 us | 34.7x | + +The Gosper iterator on flat (same-resolution) cells is only ~10-15% faster than the direct approach, since both enumerate the same set of edges. The real win comes from compaction: `compactCells` replaces complete child sets with their parents, and the Gosper iterator then walks the parent's boundary directly, skipping all internal edges. The speedup grows with resolution because compaction removes exponentially more cells. + +## Remaining work + +- Decide on public API surface for `cellsToMultiPolygonGosper` (currently internal) +- Add tests for the Gosper integration path (error cases, compacted input, mixed resolutions) +- Add input validation to `cellsToMultiPolygonGosper`: require sorted, pre-compacted input and verify cheaply with a linear scan (sorted order, no overlapping ancestors, no complete sibling sets). May need right-to-left scan for the sibling check. +- Consider whether the public `cellsToMultiPolygon` should auto-compact internally and delegate to the Gosper path diff --git a/edge_iter_plan.md b/edge_iter_plan.md new file mode 100644 index 0000000000..1932b3e0e9 --- /dev/null +++ b/edge_iter_plan.md @@ -0,0 +1,145 @@ +# GosperIter: Fresh Branch Plan + +Reimplement the EdgeIter from the `child_edge_iter` branch on a clean +branch off `master` as `GosperIter`, incorporating improvements observed +in the Zig implementation (`h3cellset/src/edge_iter.zig`). + +The name `GosperIter` reflects the geometry: the iterator traces the +Gosper island outline (the boundary of a cell's child set), distinguishing +it from a hypothetical iterator that yields all directed edges (interior ++ boundary). + +## Background + +The `child_edge_iter` branch has diverged from master and carries scratch +files (`notes.c`, `todo.md`, justfile changes). The algorithm itself is +solid -- identical to the proven Zig version -- but the C code misses a +couple of cleanups the Zig version got right. + +## Naming Convention + +Follow the existing iterator naming pattern from `iterators.h`: + +| Existing | GosperIter | +|------------------------|-------------------------------| +| `IterCellsChildren` | `IterGosper` | +| `iterInitParent(h, r)` | `iterInitGosper(h, childRes)` | +| `iterStepChild(&iter)` | `iterStepGosper(&iter)` | +| `.h` (current cell) | `.e` (current directed edge) | + +The old branch used `init_EdgeIter` / `step_EdgeIter` which doesn't +match codebase conventions. + +## Scope + +- **Internal API only** -- no `DECLSPEC`, no `H3_EXPORT`. This is a + building block for `cellsToMultiPolygon`, not a public function. +- **No input validation** -- callers are internal and already validate + upstream. +- **Own header/source** -- `gosperIter.h` / `gosperIter.c`, separate + from `iterators.h` since the concern is different. + +## Priority: Clarity Over Speed + +This iterator is expected to provide a significant speed benefit over +the current approach (enumerating child cells then computing boundaries), +but the goal of this initial implementation is **clarity for code +review**, not maximum performance. An optimization pass will follow. + +That said, avoid structural choices that would require large refactors +later. Specifically: keep the iterator state self-contained in the +struct (no heap allocation), keep the per-step work O(1), and keep the +boundary walk logic in its own static function so it can be replaced or +inlined independently. The current algorithm already satisfies all of +these -- just don't trade them away for readability shortcuts. + +## Files to Create + +| File | Description | +|---------------------------------------------|-------------------------| +| `src/h3lib/include/gosperIter.h` | Struct + init/step API | +| `src/h3lib/lib/gosperIter.c` | Iterator implementation | +| `src/apps/testapps/testGosperIter.c` | Test suite | +| `src/apps/benchmarks/benchmarkGosperIter.c` | Benchmark suite | + +## Files to Modify + +| File | Change | +|--------------------------------------|---------------------------------------| +| `CMakeLists.txt` | Add header, source, benchmark entries | +| `CMakeTests.cmake` | Add `testGosperIter` test entry | +| `src/h3lib/include/mathExtensions.h` | Add `MIN` macro | + +## Improvements Over Current Branch + +### 1. Same-resolution fast path in `init` + +Borrowed from Zig. When `parent_res == child_res`, skip the `I` array +setup, `H3_SET_RESOLUTION`, and digit-setting entirely. Just set mode +and reserved bits on the input cell directly. + +```c +if (parent_res == child_res) { + // No digit setup, resolution change, or I[] init needed + iter.e = h; + H3_SET_MODE(iter.e, H3_DIRECTEDEDGE_MODE); + H3_SET_RESERVED_BITS(iter.e, edge_seq[0]); + iter.num_edges = is_pentagon ? 5 : 6; + return iter; +} +``` + +### 2. Simpler same-resolution stepping + +Borrowed from Zig. Replace the cryptic modular-arithmetic loop: + +```c +// Before (current branch) +do { + ei->i = (ei->i + 7) % 6; +} while (ei->is_pentagon && edge_seq[ei->i] == 1); +``` + +With a straight increment and pentagon skip: + +```c +// After +iter->_i++; +if (iter->_isPentagon && iter->_i == 1) iter->_i = 2; +``` + +This works because `num_edges` guarantees we never wrap past index 5. +Clearer intent, no modular arithmetic for the simple case. + +### 3. Consistent naming + +Rename `init_EdgeIter` / `step_EdgeIter` to `iterInitGosper` / +`iterStepGosper`. Prefix internal struct fields with `_` to match +`IterCellsChildren` convention (e.g. `_parentRes`, `_skipDigit`). + +### 4. Drop scratch files + +The new branch carries only the iterator itself -- no `notes.c`, +`todo.md`, `.gitignore` changes, or `justfile` additions. + +## Not Changing + +- **API shape**: `init` returns struct, `.e` is current value (`H3_NULL` + when done), `step` advances. Matches existing `IterCellsChildren` + convention. +- **Core algorithm**: The boundary walk (`base[18]`, `edge_seq[6]`, + `+19 mod 18`, `-6` parent adjustment, recursive `step_boundaryCell`) + is proven correct in both C and Zig with exhaustive testing. +- **Test coverage**: Exhaustive res 0-4 + specific edge sequences + + loop connectivity + per-edge validity. Already matches the Zig suite. +- **Benchmark coverage**: Same-res through +11 gap, hex and pentagon, + multiple base resolutions. + +## Steps + +1. `git checkout master && git checkout -b gosper_iter` +2. Create `gosperIter.h` and `gosperIter.c` with improvements above +3. Create `testGosperIter.c` and `benchmarkGosperIter.c` +4. Update `CMakeLists.txt` and `CMakeTests.cmake` +5. Add `MIN` macro to `mathExtensions.h` +6. Build and run tests: `just build && just test` diff --git a/justfile b/justfile new file mode 100644 index 0000000000..048513fbb4 --- /dev/null +++ b/justfile @@ -0,0 +1,64 @@ +init: purge + mkdir build + +build: + cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + +purge: + rm -rf build + rm -rf *.trace + rm -rf .ipynb_checkpoints + rm -rf .cache + rm -rf .claude + +test-fast: build + cd build; make test-fast + +test-slow: build + cd build; make test + +# Run a single test binary. Dots (progress) go to /dev/null; failures print to stderr. +test-one TEST: build + ./build/bin/{{TEST}} > /dev/null + +test: + just test-fast + +bench: build + ./build/bin/benchmarkCellsToPolyAlgos + +coverage: purge + mkdir build + cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage + open build/coverage/index.html + +# Show uncovered lines/branches for a source file from the lcov .info data. +# Run `just coverage` first. Example: just coverage-gaps linkedGeo.c +coverage-gaps FILE: + #!/usr/bin/env bash + set -e + info="build/coverage.cleaned.info" + if [ ! -f "$info" ]; then + echo "No coverage data found — run 'just coverage' first." + exit 1 + fi + # Extract the section for this file + section=$(sed -n "/SF:.*\/{{FILE}}$/,/end_of_record/p" "$info") + if [ -z "$section" ]; then + echo "File {{FILE}} not found in coverage data." + echo "Available files:" + grep '^SF:' "$info" | sed 's|.*src/h3lib/||' + exit 1 + fi + echo "=== Summary ===" + lf=$(echo "$section" | grep '^LF:' | cut -d: -f2) + lh=$(echo "$section" | grep '^LH:' | cut -d: -f2) + brf=$(echo "$section" | grep '^BRF:' | cut -d: -f2) + brh=$(echo "$section" | grep '^BRH:' | cut -d: -f2) + echo "Lines: $lh/$lf Branches: $brh/$brf" + echo "" + echo "=== Uncovered lines (DA:line,0) ===" + echo "$section" | grep '^DA:' | awk -F'[,:]' '$3 == 0 {print " line " $2}' || echo " (none)" + echo "" + echo "=== Untaken branches (BRDA:line,block,branch,0) ===" + echo "$section" | grep '^BRDA:' | awk -F'[,:]' '$5 == 0 {print " line " $2 " branch " $4}' || echo " (none)" diff --git a/src/apps/applib/include/benchmark.h b/src/apps/applib/include/benchmark.h index 276b1cd4c2..825808e292 100644 --- a/src/apps/applib/include/benchmark.h +++ b/src/apps/applib/include/benchmark.h @@ -30,6 +30,7 @@ #ifdef _WIN32 #define NOGDI +#define WIN32_LEAN_AND_MEAN #include #define START_TIMER \ diff --git a/src/apps/benchmarks/benchmarkCellsToPolyAlgos.c b/src/apps/benchmarks/benchmarkCellsToPolyAlgos.c index e6a079d597..38f26a7fed 100644 --- a/src/apps/benchmarks/benchmarkCellsToPolyAlgos.c +++ b/src/apps/benchmarks/benchmarkCellsToPolyAlgos.c @@ -41,6 +41,76 @@ }); \ } while (0) +#define BENCHMARK_GOSPER(NAME, ITERS, RES) \ + do { \ + BENCHMARK(gosper_##NAME, ITERS, { \ + GeoMultiPolygon mpoly; \ + cellsToMultiPolygonGosper(cells, numCells, RES, &mpoly); \ + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); \ + }); \ + } while (0) + +#define BENCHMARK_GOSPER_COMPACT(NAME, ITERS, CELLS, N, RES) \ + do { \ + BENCHMARK(gosper_compact_##NAME, ITERS, { \ + GeoMultiPolygon mpoly; \ + cellsToMultiPolygonGosper(CELLS, N, RES, &mpoly); \ + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); \ + }); \ + } while (0) + +// Helper: fill Colorado cells at a given resolution. +// Caller must free *outCells. +void coloradoCells(int res, H3Index **outCells, int64_t *outNumCells) { + LatLng verts[] = { + {37.0, -109.0}, + {37.0, -102.0}, + {41.0, -102.0}, + {41.0, -109.0}, + }; + for (int i = 0; i < 4; i++) { + verts[i].lat = H3_EXPORT(degsToRads)(verts[i].lat); + verts[i].lng = H3_EXPORT(degsToRads)(verts[i].lng); + } + GeoPolygon polygon = { + .numHoles = 0, + .holes = NULL, + .geoloop = {.numVerts = 4, .verts = verts}, + }; + + int64_t maxCells; + H3_EXPORT(maxPolygonToCellsSize)(&polygon, res, 0, &maxCells); + + H3Index *cells = calloc(maxCells, sizeof(H3Index)); + H3_EXPORT(polygonToCells)(&polygon, res, 0, cells); + + int64_t numCells = 0; + for (int64_t i = 0; i < maxCells; i++) { + if (cells[i] != H3_NULL) { + cells[numCells++] = cells[i]; + } + } + cells = realloc(cells, numCells * sizeof(H3Index)); + + *outCells = cells; + *outNumCells = numCells; +} + +// Helper: compact a cell set. Caller must free *outCompacted. +void compactAndCount(H3Index *cells, int64_t numCells, H3Index **outCompacted, + int64_t *outNumCompacted) { + H3Index *compacted = calloc(numCells, sizeof(H3Index)); + H3_EXPORT(compactCells)(cells, compacted, numCells); + int64_t numCompacted = 0; + for (int64_t i = 0; i < numCells; i++) { + if (compacted[i] != H3_NULL) { + compacted[numCompacted++] = compacted[i]; + } + } + *outCompacted = compacted; + *outNumCompacted = numCompacted; +} + BEGIN_BENCHMARKS(); { @@ -108,45 +178,54 @@ BEGIN_BENCHMARKS(); } { - // Square approximating Colorado. - // (4 corners, counterclockwise from southwest) int res = 6; - LatLng verts[] = { - {37.0, -109.0}, - {37.0, -102.0}, - {41.0, -102.0}, - {41.0, -109.0}, - }; + H3Index *cells; + int64_t numCells; + coloradoCells(res, &cells, &numCells); - for (int i = 0; i < 4; i++) { - verts[i].lat = H3_EXPORT(degsToRads)(verts[i].lat); - verts[i].lng = H3_EXPORT(degsToRads)(verts[i].lng); - } + BENCHMARK_LINKED(colorado, 100); + BENCHMARK_DIRECT(colorado, 100); + BENCHMARK_GOSPER(colorado, 100, res); - GeoPolygon polygon = { - .numHoles = 0, - .holes = NULL, - .geoloop = {.numVerts = 4, .verts = verts}, - }; + H3Index *compacted; + int64_t numCompacted; + compactAndCount(cells, numCells, &compacted, &numCompacted); + BENCHMARK_GOSPER_COMPACT(colorado, 100, compacted, numCompacted, res); + free(compacted); + free(cells); +} - int64_t maxCells; - H3_EXPORT(maxPolygonToCellsSize)(&polygon, res, 0, &maxCells); +{ + int res = 7; + H3Index *cells; + int64_t numCells; + coloradoCells(res, &cells, &numCells); - H3Index *cells = calloc(maxCells, sizeof(H3Index)); - H3_EXPORT(polygonToCells)(&polygon, res, 0, cells); + BENCHMARK_DIRECT(colorado7, 10); + BENCHMARK_GOSPER(colorado7, 10, res); - // Move all valid cells to the front and resize - int64_t numCells = 0; - for (int64_t i = 0; i < maxCells; i++) { - if (cells[i] != H3_NULL) { - cells[numCells++] = cells[i]; - } - } - cells = realloc(cells, numCells * sizeof(H3Index)); + H3Index *compacted; + int64_t numCompacted; + compactAndCount(cells, numCells, &compacted, &numCompacted); + BENCHMARK_GOSPER_COMPACT(colorado7, 10, compacted, numCompacted, res); + free(compacted); + free(cells); +} - BENCHMARK_LINKED(colorado, 100); - BENCHMARK_DIRECT(colorado, 100); +{ + int res = 8; + H3Index *cells; + int64_t numCells; + coloradoCells(res, &cells, &numCells); + + BENCHMARK_DIRECT(colorado8, 1); + BENCHMARK_GOSPER(colorado8, 1, res); + H3Index *compacted; + int64_t numCompacted; + compactAndCount(cells, numCells, &compacted, &numCompacted); + BENCHMARK_GOSPER_COMPACT(colorado8, 1, compacted, numCompacted, res); + free(compacted); free(cells); } diff --git a/src/apps/benchmarks/benchmarkGosperIter.c b/src/apps/benchmarks/benchmarkGosperIter.c new file mode 100644 index 0000000000..98232b708a --- /dev/null +++ b/src/apps/benchmarks/benchmarkGosperIter.c @@ -0,0 +1,53 @@ +/* + * Copyright 2026 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "benchmark.h" +#include "iterators.h" + +H3Index hex2 = 0x820887fffffffff; // res 2 hexagon +H3Index pent2 = 0x820807fffffffff; // res 2 pentagon + +void exhaustGosper(H3Index cell, int childRes) { + IterEdgesGosper iter = iterInitGosper(cell, childRes); + while (iter.e) iterStepGosper(&iter); +} + +BEGIN_BENCHMARKS(); + +// res: 2 -> 2 +BENCHMARK(hex_plus0, 50000000, { exhaustGosper(hex2, 2); }); +BENCHMARK(pent_plus0, 50000000, { exhaustGosper(pent2, 2); }); + +// res: 2 -> 3 +BENCHMARK(hex_plus1, 500000, { exhaustGosper(hex2, 3); }); +BENCHMARK(pent_plus1, 500000, { exhaustGosper(pent2, 3); }); + +// res: 2 -> 4 +BENCHMARK(hex_plus2, 100000, { exhaustGosper(hex2, 4); }); +BENCHMARK(pent_plus2, 100000, { exhaustGosper(pent2, 4); }); + +// res: 2 -> 7 +BENCHMARK(hex_plus5, 1000, { exhaustGosper(hex2, 7); }); +BENCHMARK(pent_plus5, 1000, { exhaustGosper(pent2, 7); }); + +// res: 2 -> 10 +BENCHMARK(hex_plus8, 100, { exhaustGosper(hex2, 10); }); +BENCHMARK(pent_plus8, 100, { exhaustGosper(pent2, 10); }); + +// res: 2 -> 13 +BENCHMARK(hex_plus11, 10, { exhaustGosper(hex2, 13); }); +BENCHMARK(pent_plus11, 10, { exhaustGosper(pent2, 13); }); + +END_BENCHMARKS(); diff --git a/src/apps/testapps/testGosperIter.c b/src/apps/testapps/testGosperIter.c new file mode 100644 index 0000000000..59a09d8968 --- /dev/null +++ b/src/apps/testapps/testGosperIter.c @@ -0,0 +1,320 @@ +/* + * Copyright 2026 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include + +#include "iterators.h" +#include "mathExtensions.h" +#include "test.h" +#include "utility.h" + +// Check iterator mechanics: +// 1. produces correct number of edges +// 2. each step produces an edge that's different from the previous edge +// 3. stays H3_NULL after finishing +void check_iterator_mechanics(H3Index h, int childRes) { + bool isPent = H3_EXPORT(isPentagon)(h); + int parentRes = H3_EXPORT(getResolution)(h); + int numSides = isPent ? 5 : 6; + int64_t expectedEdges = numSides * _ipow(3, childRes - parentRes); + + IterEdgesGosper iter = iterInitGosper(h, childRes); + t_assert(iter.remaining == expectedEdges, "correct number of edges"); + + int64_t i; + for (i = 0; iter.e; i++) { + H3Index prev = iter.e; + iterStepGosper(&iter); + t_assert(iter.e != prev, "edge should advance"); + } + t_assert(i == expectedEdges, "correct number of edges"); + + // iterator should stay exhausted + for (int j = 0; j < 100; j++) { + t_assert(iter.remaining == 0, "num edges to zero"); + t_assert(iter.e == H3_NULL, "iterator exhausted"); + iterStepGosper(&iter); + } +} + +// Per-edge correctness: +// 1. each edge is a valid directed edge +// 2. each edge is at childRes +// 3. confirm the edge is on the Gosper island boundary: +// origin cell is a child of parent cell, destination cell is not +void check_edge_correctness(H3Index h, int childRes) { + int parentRes = H3_EXPORT(getResolution)(h); + IterEdgesGosper iter = iterInitGosper(h, childRes); + + while (iter.e) { + t_assert(H3_EXPORT(isValidDirectedEdge)(iter.e), "edge is valid"); + t_assert(H3_EXPORT(getResolution)(iter.e) == childRes, + "correct resolution"); + + H3Index s; // source/origin cell + H3Index d; // destination cell + H3Index ps; // parent of source + H3Index pd; // parent of destination + + t_assertSuccess(H3_EXPORT(getDirectedEdgeOrigin)(iter.e, &s)); + t_assertSuccess(H3_EXPORT(getDirectedEdgeDestination)(iter.e, &d)); + t_assertSuccess(H3_EXPORT(cellToParent)(s, parentRes, &ps)); + t_assertSuccess(H3_EXPORT(cellToParent)(d, parentRes, &pd)); + + t_assert(ps == h, "origin *is* child of parent cell"); + t_assert(pd != h, "destination *is not* child of parent cell"); + + iterStepGosper(&iter); + } +} + +// Check that edge_a's last vertex matches edge_b's first vertex, +// i.e. the edges connect end-to-start. Tolerance is relative to +// the shorter edge length to handle varying resolutions. +// TODO: Replace floating-point distance check with exact vertex comparison +// once an edgeToVertexes function exists. +bool do_edges_connect(H3Index edge_a, H3Index edge_b) { + double len_a, len_b; + H3_EXPORT(edgeLengthRads)(edge_a, &len_a); + H3_EXPORT(edgeLengthRads)(edge_b, &len_b); + double tol = MIN(len_a, len_b) / 1000.0; + + CellBoundary bd_a, bd_b; + t_assertSuccess(H3_EXPORT(directedEdgeToBoundary)(edge_a, &bd_a)); + t_assertSuccess(H3_EXPORT(directedEdgeToBoundary)(edge_b, &bd_b)); + + LatLng end_a = bd_a.verts[bd_a.numVerts - 1]; + LatLng start_b = bd_b.verts[0]; + double dist = H3_EXPORT(greatCircleDistanceRads)(&end_a, &start_b); + + return dist < tol; +} + +// Loop structure: +// 1. consecutive edges share an endpoint (floating-point test) +// 2. loop closes (returns to where it started) +void check_loop_valid(H3Index h, int childRes) { + IterEdgesGosper iter = iterInitGosper(h, childRes); + H3Index first_edge = iter.e; + H3Index prev_edge = iter.e; + iterStepGosper(&iter); + + while (iter.e) { + t_assert(do_edges_connect(prev_edge, iter.e), "edges should connect"); + + prev_edge = iter.e; + iterStepGosper(&iter); + } + + t_assert(do_edges_connect(prev_edge, first_edge), "loop should close"); +} + +void check_gosper_island(H3Index h, int childRes) { + check_iterator_mechanics(h, childRes); + check_edge_correctness(h, childRes); + check_loop_valid(h, childRes); +} + +// Verify iterator produces some cyclic rotation of the expected edges +void check_expected_edges(H3Index h, int childRes, H3Index *expected, + int numExpected) { + IterEdgesGosper iter = iterInitGosper(h, childRes); + t_assert(iter.remaining == numExpected, "correct number of edges"); + + // Find rotation offset: where does the iterator's first edge appear? + int offset = -1; + for (int i = 0; i < numExpected; i++) { + if (expected[i] == iter.e) { + offset = i; + break; + } + } + t_assert(offset >= 0, "first edge found in expected array"); + + for (int i = 0; iter.e; i++) { + t_assert(iter.e == expected[(offset + i) % numExpected], + "correct edge"); + iterStepGosper(&iter); + } +} + +SUITE(gosper_iter) { + // Exhaustive test of all cells, resolutions 0,1,2,3,4. + // Expanding to Gosper islands up to res 4. + TEST(coarse_resolutions_exhaustive) { + for (int childRes = 0; childRes <= 4; childRes++) { + for (int parentRes = 0; parentRes <= childRes; parentRes++) { + IterCellsResolution cell_iter = iterInitRes(parentRes); + while (cell_iter.h) { + check_gosper_island(cell_iter.h, childRes); + iterStepRes(&cell_iter); + } + } + } + } + + // Since an exhaustive test of all resolutions would take too long, + // we test a few specific cells at finer resolutions, including + // expanding up to res 15. + TEST(finer_resolutions) { + H3Index cells[] = { + 0x8508000ffffffff, // res 5 hexagon + 0x85080003fffffff, // res 5 pentagon + + 0x8808000009fffff, // res 8 hexagon + 0x8808000001fffff, // res 8 pentagon + + 0x8a0800000017fff, // res 10 hexagon + 0x8a0800000007fff, // res 10 pentagon + + 0x8e754e64992d6c7, // res 14 hexagon + 0x8e0800000000007, // res 14 pentagon + + 0x8f754e64992d6d8, // res 15 hexagon + 0x8f0800000000000, // res 15 pentagon + }; + for (int i = 0; i < (int)ARRAY_SIZE(cells); i++) { + int parentRes = H3_EXPORT(getResolution)(cells[i]); + for (int gap = 0; gap <= 5; gap++) { + int childRes = parentRes + gap; + if (childRes > 15) continue; + check_gosper_island(cells[i], childRes); + } + } + } + + // `check_expected_edges` tests below check for a specific sequence + // edges for a given cell and resolution input. Any cyclic rotation + // should pass the test. + + TEST(just_hex_2) { + H3Index h = 0x820887fffffffff; + int res = 2; + H3Index expected[] = { + 0x1320887fffffffff, 0x1120887fffffffff, 0x1520887fffffffff, + 0x1420887fffffffff, 0x1620887fffffffff, 0x1220887fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(just_pent_2) { + H3Index h = 0x820807fffffffff; + int res = 2; + H3Index expected[] = { + 0x1320807fffffffff, 0x1520807fffffffff, 0x1420807fffffffff, + 0x1620807fffffffff, 0x1220807fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(pent_edges_2_to_4) { + H3Index h = 0x820807fffffffff; + int res = 4; + H3Index expected[] = { + 0x11408053ffffffff, 0x15408053ffffffff, 0x1140805bffffffff, + 0x1540805bffffffff, 0x1440805bffffffff, 0x15408059ffffffff, + 0x14408059ffffffff, 0x16408059ffffffff, 0x1440805dffffffff, + 0x1540804bffffffff, 0x1440804bffffffff, 0x15408049ffffffff, + 0x14408049ffffffff, 0x16408049ffffffff, 0x1440804dffffffff, + 0x1640804dffffffff, 0x1240804dffffffff, 0x16408045ffffffff, + 0x14408069ffffffff, 0x16408069ffffffff, 0x1440806dffffffff, + 0x1640806dffffffff, 0x1240806dffffffff, 0x16408065ffffffff, + 0x12408065ffffffff, 0x13408065ffffffff, 0x12408067ffffffff, + 0x1640802dffffffff, 0x1240802dffffffff, 0x16408025ffffffff, + 0x12408025ffffffff, 0x13408025ffffffff, 0x12408027ffffffff, + 0x13408027ffffffff, 0x11408027ffffffff, 0x13408023ffffffff, + 0x12408035ffffffff, 0x13408035ffffffff, 0x12408037ffffffff, + 0x13408037ffffffff, 0x11408037ffffffff, 0x13408033ffffffff, + 0x11408033ffffffff, 0x15408033ffffffff, 0x1140803bffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(pent_edges_1_to_3) { + H3Index h = 0x81083ffffffffff; + int res = 3; + H3Index expected[] = { + 0x113082bfffffffff, 0x1330829fffffffff, 0x1130829fffffffff, + 0x1530829fffffffff, 0x113082dfffffffff, 0x153082dfffffffff, + 0x143082dfffffffff, 0x153082cfffffffff, 0x143082cfffffffff, + 0x1530821fffffffff, 0x1130825fffffffff, 0x1530825fffffffff, + 0x1430825fffffffff, 0x1530824fffffffff, 0x1430824fffffffff, + 0x1630824fffffffff, 0x1430826fffffffff, 0x1630826fffffffff, + 0x1430835fffffffff, 0x1530834fffffffff, 0x1430834fffffffff, + 0x1630834fffffffff, 0x1430836fffffffff, 0x1630836fffffffff, + 0x1230836fffffffff, 0x1630832fffffffff, 0x1230832fffffffff, + 0x1630814fffffffff, 0x1430816fffffffff, 0x1630816fffffffff, + 0x1230816fffffffff, 0x1630812fffffffff, 0x1230812fffffffff, + 0x1330812fffffffff, 0x1230813fffffffff, 0x1330813fffffffff, + 0x123081efffffffff, 0x163081afffffffff, 0x123081afffffffff, + 0x133081afffffffff, 0x123081bfffffffff, 0x133081bfffffffff, + 0x113081bfffffffff, 0x1330819fffffffff, 0x1130819fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(hex_edges_1_to_3) { + H3Index h = 0x81c67ffffffffff; + int res = 3; + H3Index expected[] = { + 0x133c64afffffffff, 0x123c64bfffffffff, 0x133c64bfffffffff, + 0x113c64bfffffffff, 0x133c649fffffffff, 0x113c649fffffffff, + 0x153c649fffffffff, 0x113c64dfffffffff, 0x153c64dfffffffff, + 0x113c66bfffffffff, 0x133c669fffffffff, 0x113c669fffffffff, + 0x153c669fffffffff, 0x113c66dfffffffff, 0x153c66dfffffffff, + 0x143c66dfffffffff, 0x153c66cfffffffff, 0x143c66cfffffffff, + 0x153c661fffffffff, 0x113c665fffffffff, 0x153c665fffffffff, + 0x143c665fffffffff, 0x153c664fffffffff, 0x143c664fffffffff, + 0x163c664fffffffff, 0x143c666fffffffff, 0x163c666fffffffff, + 0x143c675fffffffff, 0x153c674fffffffff, 0x143c674fffffffff, + 0x163c674fffffffff, 0x143c676fffffffff, 0x163c676fffffffff, + 0x123c676fffffffff, 0x163c672fffffffff, 0x123c672fffffffff, + 0x163c654fffffffff, 0x143c656fffffffff, 0x163c656fffffffff, + 0x123c656fffffffff, 0x163c652fffffffff, 0x123c652fffffffff, + 0x133c652fffffffff, 0x123c653fffffffff, 0x133c653fffffffff, + 0x123c65efffffffff, 0x163c65afffffffff, 0x123c65afffffffff, + 0x133c65afffffffff, 0x123c65bfffffffff, 0x133c65bfffffffff, + 0x113c65bfffffffff, 0x133c659fffffffff, 0x113c659fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(hex_edges_2_to_4) { + H3Index h = 0x82c64ffffffffff; + int res = 4; + H3Index expected[] = { + 0x134c6497ffffffff, 0x114c6497ffffffff, 0x134c6493ffffffff, + 0x114c6493ffffffff, 0x154c6493ffffffff, 0x114c649bffffffff, + 0x154c649bffffffff, 0x144c649bffffffff, 0x154c6499ffffffff, + 0x114c64d3ffffffff, 0x154c64d3ffffffff, 0x114c64dbffffffff, + 0x154c64dbffffffff, 0x144c64dbffffffff, 0x154c64d9ffffffff, + 0x144c64d9ffffffff, 0x164c64d9ffffffff, 0x144c64ddffffffff, + 0x154c64cbffffffff, 0x144c64cbffffffff, 0x154c64c9ffffffff, + 0x144c64c9ffffffff, 0x164c64c9ffffffff, 0x144c64cdffffffff, + 0x164c64cdffffffff, 0x124c64cdffffffff, 0x164c64c5ffffffff, + 0x144c64e9ffffffff, 0x164c64e9ffffffff, 0x144c64edffffffff, + 0x164c64edffffffff, 0x124c64edffffffff, 0x164c64e5ffffffff, + 0x124c64e5ffffffff, 0x134c64e5ffffffff, 0x124c64e7ffffffff, + 0x164c64adffffffff, 0x124c64adffffffff, 0x164c64a5ffffffff, + 0x124c64a5ffffffff, 0x134c64a5ffffffff, 0x124c64a7ffffffff, + 0x134c64a7ffffffff, 0x114c64a7ffffffff, 0x134c64a3ffffffff, + 0x124c64b5ffffffff, 0x134c64b5ffffffff, 0x124c64b7ffffffff, + 0x134c64b7ffffffff, 0x114c64b7ffffffff, 0x134c64b3ffffffff, + 0x114c64b3ffffffff, 0x154c64b3ffffffff, 0x114c64bbffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } +} diff --git a/src/h3lib/include/cellsToMultiPoly.h b/src/h3lib/include/cellsToMultiPoly.h index f22af2dc8f..4330eded13 100644 --- a/src/h3lib/include/cellsToMultiPoly.h +++ b/src/h3lib/include/cellsToMultiPoly.h @@ -220,6 +220,27 @@ static inline void destroySortablePolyVerts(SortablePoly *spolys, } } +/** + * Create a GeoMultiPolygon from a set of (possibly compacted) H3 cells + * using the Gosper island boundary iterator. Much faster than the flat + * cellsToMultiPolygon when input has been run through compactCells first. + * + * Unlike cellsToMultiPolygon, this function does NOT validate, deduplicate, + * or compact the input. The caller must ensure cells are valid, non- + * overlapping, and at resolutions <= targetRes. For best performance, + * compact the input with compactCells before calling. + * + * @param cells Array of valid H3 cells (may be at mixed resolutions). + * @param numCells Number of cells. + * @param targetRes Resolution of the output polygon edges. Must be >= + * the resolution of every cell in the array. + * @param out Output GeoMultiPolygon. Caller frees with + * destroyGeoMultiPolygon. + * @return E_SUCCESS on success + */ +H3Error cellsToMultiPolygonGosper(const H3Index *cells, int64_t numCells, + int targetRes, GeoMultiPolygon *out); + /** @brief Create a GeoMultiPolygon from a set of cells * * NOTE: This definition is tentative as we work to finish the implementation. diff --git a/src/h3lib/include/iterators.h b/src/h3lib/include/iterators.h index 9e1de9bfdc..899e08c9eb 100644 --- a/src/h3lib/include/iterators.h +++ b/src/h3lib/include/iterators.h @@ -22,8 +22,10 @@ #ifndef ITERATORS_H #define ITERATORS_H +#include #include +#include "h3Index.h" #include "h3api.h" /** @@ -84,4 +86,45 @@ typedef struct { DECLSPEC IterCellsResolution iterInitRes(int res); DECLSPEC void iterStepRes(IterCellsResolution *iter); +/** + * IterEdgesGosper: iterator for the directed edges on the boundary of a cell's + * child set (the Gosper island outline) at a given resolution. + * + * Each yielded edge is a directed edge whose origin cell is a child of the + * parent cell and whose destination cell is not. The edges form a closed loop + * in counter-clockwise order around the Gosper island. The starting edge is + * not guaranteed; any cyclic rotation of the sequence is a valid output. + * + * Constructor: + * + * Initialize with `iterInitGosper(cell, childRes)`. + * Requires `childRes >= getResolution(cell)` and a valid cell. + * + * Iteration: + * + * Step iterator with `iterStepGosper`. + * The current edge is accessed via the `.e` member. + * The `.remaining` member gives the number of edges left to yield + * (including the current edge). + * When the iterator is exhausted, `.e` will be `H3_NULL` even after + * further calls to `iterStepGosper`. + */ +typedef struct { + H3Index e; // Current directed edge (H3_NULL when exhausted) + int64_t remaining; // Number of edges left to yield (including current) + + // Internal state — cyclic position trackers: + int8_t _edgePos; // Position in edge_dir cycle (period 6) + int8_t _walkPos[MAX_H3_RES + 1]; // Position in walk_digit cycle + // (period 18) per resolution level + + // Additional internal state + int8_t _parentRes; + int8_t _childRes; + bool _isPentagon; +} IterEdgesGosper; + +DECLSPEC IterEdgesGosper iterInitGosper(H3Index h, int childRes); +DECLSPEC void iterStepGosper(IterEdgesGosper *iter); + #endif diff --git a/src/h3lib/include/mathExtensions.h b/src/h3lib/include/mathExtensions.h index 1ac6cfdba3..c5715b2c38 100644 --- a/src/h3lib/include/mathExtensions.h +++ b/src/h3lib/include/mathExtensions.h @@ -28,6 +28,11 @@ */ #define MAX(a, b) (((a) > (b)) ? (a) : (b)) +/** + * MIN returns the minimum of two values. + */ +#define MIN(a, b) (((a) < (b)) ? (a) : (b)) + /** Evaluates to true if a + b would overflow for int32 */ static inline bool ADD_INT32S_OVERFLOWS(int32_t a, int32_t b) { if (a > 0) { diff --git a/src/h3lib/lib/cellsToMultiPoly.c b/src/h3lib/lib/cellsToMultiPoly.c index 3145cee5d5..d9aafb4a2d 100644 --- a/src/h3lib/lib/cellsToMultiPoly.c +++ b/src/h3lib/lib/cellsToMultiPoly.c @@ -10,6 +10,8 @@ #include "area.h" #include "h3Assert.h" #include "h3api.h" +#include "iterators.h" +#include "mathExtensions.h" static inline H3Error validateCellSet(const H3Index *cells, const int64_t numCells) { @@ -572,6 +574,38 @@ static H3Error createMultiPolygon(SortableLoopSet loopset, return E_SUCCESS; } +/** + * Shared pipeline: cancel arc pairs, extract loops, build polygons. + * Takes ownership of arcset on success (caller must not free). + * On error, arcset is cleaned up here. + */ +static H3Error arcSetToMultiPolygon(ArcSet *arcset, GeoMultiPolygon *out) { + H3Error err = cancelArcPairs(*arcset); + if (NEVER(err)) { + destroyArcSet(arcset); + return err; + } + + SortableLoopSet loopset; + err = createSortableLoopSet(*arcset, &loopset); + if (err) { + destroyArcSet(arcset); + return err; + } + + err = createMultiPolygon(loopset, out); + if (err) { + destroySortableLoopSet(&loopset); + destroyArcSet(arcset); + return err; + } + + destroyArcSet(arcset); + destroySortableLoopSetShallow(&loopset); + + return E_SUCCESS; +} + /** * Create a GeoMultiPolygon from a set of H3 cells. * @@ -619,47 +653,142 @@ H3Error H3_EXPORT(cellsToMultiPolygon)(const H3Index *cells, return E_SUCCESS; } - // arcset initializes with separate doubly-linked loops for each cell, - // each in their own connected component ArcSet arcset; err = createArcSet(cells, numCells, &arcset); if (err) return err; - // Cancel out pairs of edges, updating the doubly-linked loops and merging - // them into a single connected component - err = cancelArcPairs(arcset); - if (NEVER(err)) { - destroyArcSet(&arcset); - return err; + return arcSetToMultiPolygon(&arcset, out); +} + +/** + * Count total Gosper boundary edges across a set of (possibly compacted) + * cells, each expanded to targetRes. + */ +static int64_t countGosperEdges(const H3Index *cells, int64_t numCells, + int targetRes) { + int64_t total = 0; + for (int64_t i = 0; i < numCells; i++) { + int parentRes = H3_GET_RESOLUTION(cells[i]); + bool pent = H3_EXPORT(isPentagon)(cells[i]); + int faces = pent ? 5 : 6; + total += faces * _ipow(3, targetRes - parentRes); } + return total; +} - /* - Extract all loops and sort them by: - 1) their connected component, and then by - 2) the loop area. - This makes loops for each polygon contiguous in memory. - Within each polygon, the sorting makes the loop with the smallest enclosed - area come first (accounting for loop winding direction), - which is what we take to be the outer loop for that polygon. - */ - SortableLoopSet loopset; - err = createSortableLoopSet(arcset, &loopset); - if (err) { - destroyArcSet(&arcset); - return err; +/** + * Create an ArcSet from Gosper boundary edges of (possibly compacted) cells. + * + * For each cell, iterInitGosper yields edges in geometric (CCW) order around + * that cell's Gosper island at targetRes. Edges from different cells that + * share a boundary will be paired and canceled later by cancelArcPairs. + */ +static H3Error createArcSetGosper(const H3Index *cells, int64_t numCells, + int targetRes, ArcSet *arcset) { + int64_t numArcs = countGosperEdges(cells, numCells, targetRes); + int64_t numBuckets = numArcs * HASH_TABLE_MULTIPLIER; + + arcset->numArcs = numArcs; + arcset->numBuckets = numBuckets; + arcset->arcs = H3_MEMORY(malloc)(numArcs * sizeof(Arc)); + if (!arcset->arcs) { + return E_MEMORY_ALLOC; } - // Extract polygons, since loops are contiguous in SortableLoopSet memory. - // Polygons sorted by outer loop area, decreasing. - err = createMultiPolygon(loopset, out); - if (err) { - destroySortableLoopSet(&loopset); - destroyArcSet(&arcset); - return err; + arcset->buckets = H3_MEMORY(calloc)(numBuckets, sizeof(Arc *)); + if (!arcset->buckets) { + destroyArcSet(arcset); + return E_MEMORY_ALLOC; } - destroyArcSet(&arcset); - destroySortableLoopSetShallow(&loopset); + // Fill arcs from Gosper boundary edges. Each cell's edges form a + // doubly-linked loop in their own connected component. + int64_t j = 0; + for (int64_t i = 0; i < numCells; i++) { + IterEdgesGosper iter = iterInitGosper(cells[i], targetRes); + int64_t loopStart = j; + + while (iter.e) { + arcset->arcs[j].id = iter.e; + arcset->arcs[j].isRemoved = false; + arcset->arcs[j].isVisited = false; + + // Union-find: all arcs in this cell's loop share first arc as root + arcset->arcs[j].parent = &arcset->arcs[loopStart]; + arcset->arcs[j].rank = 1; + + // Prev/next: iterator yields edges in CCW order, so link + // sequentially. We'll close the loop after the while loop. + if (j > loopStart) { + arcset->arcs[j].prev = &arcset->arcs[j - 1]; + arcset->arcs[j - 1].next = &arcset->arcs[j]; + } + + j++; + iterStepGosper(&iter); + } + + // Close the doubly-linked loop + if (j > loopStart) { + arcset->arcs[loopStart].prev = &arcset->arcs[j - 1]; + arcset->arcs[j - 1].next = &arcset->arcs[loopStart]; + } + } + + // Build hash table + for (int64_t i = 0; i < arcset->numArcs; i++) { + int64_t b = hashEdge(arcset->arcs[i].id, arcset->numBuckets); + while (arcset->buckets[b] != NULL) { + b = (b + 1) % arcset->numBuckets; + } + arcset->buckets[b] = &arcset->arcs[i]; + } return E_SUCCESS; } + +/** + * Create a GeoMultiPolygon from a set of (possibly compacted) H3 cells. + * + * Each cell is expanded to targetRes using the Gosper iterator, which yields + * only boundary edges. Cross-cell shared edges are canceled via hash lookup. + * This is faster than the flat approach for compacted input because it avoids + * enumerating and canceling the vast majority of internal edges. + * + * Unlike cellsToMultiPolygon, this function does NOT validate, deduplicate, + * or compact the input. The caller must ensure cells are valid, non- + * overlapping, and at resolutions <= targetRes. For best performance, + * compact the input with compactCells before calling. + * + * TODO: Consider requiring sorted, pre-compacted input and validating that + * cheaply. A linear scan over sorted input can verify: (1) sorted order, + * (2) no cell is an ancestor of the next (no overlaps), and (3) no + * complete set of same-resolution siblings is present (i.e., the input + * is fully compacted). Check (3) may need a right-to-left scan so that + * consecutive siblings can be counted efficiently. + * + * @param cells Array of valid H3 cells (may be at mixed resolutions). + * @param numCells Number of cells. + * @param targetRes Resolution of the output edges. Must be >= resolution + * of every cell in the array. + * @param out Output GeoMultiPolygon. Caller frees with destroyGeoMultiPolygon. + * @return E_SUCCESS on success + */ +H3Error cellsToMultiPolygonGosper(const H3Index *cells, int64_t numCells, + int targetRes, GeoMultiPolygon *out) { + if (numCells == 0) { + out->numPolygons = 0; + out->polygons = NULL; + return E_SUCCESS; + } + + H3Error err = + checkCellsToMultiPolyOverflow(numCells, HASH_TABLE_MULTIPLIER); + if (err) return err; + + ArcSet arcset; + err = createArcSetGosper(cells, numCells, targetRes, &arcset); + if (err) return err; + + return arcSetToMultiPolygon(&arcset, out); +} diff --git a/src/h3lib/lib/iterGosper.c b/src/h3lib/lib/iterGosper.c new file mode 100644 index 0000000000..b78b91e2ce --- /dev/null +++ b/src/h3lib/lib/iterGosper.c @@ -0,0 +1,231 @@ +/* + * Copyright 2026 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file iterGosper.c + * @brief Iterator for the directed edges forming the "Gosper island" outline + * of a cell's child set at a given resolution. + * + * For a cell with resolution `r`, the boundary of the cell's children at + * resolution `R` forms a closed loop of directed edges (the Gosper island). + * For a hexagon this is 6 * 3^(R - r) edges; for a pentagon, 5 * 3^(R - r). + * We refer to `R - r` as the resolution delta. + * + * ## Geometric intuition + * + * At delta 0, we just return the cell's own edges in counter-clockwise + * order, cycling through the `edge_dir` array. + * + * At delta 1 (for a hexagon), each side of the parent subdivides into 3 edges + * at the finer resolution, giving 6 * 3 = 18 edges total. The origin cells + * share the parent's digit path, extended by one digit into the child + * resolution. The `walk_digit` array maps these WALK_CYCLE_LEN walk positions + * to H3 digit values: the 6 non-center child digits in counter-clockwise order, + * each repeated 3x for the subdivision. + * + * ## Pentagon handling + * + * Pentagons have 5 sides instead of 6: the K-axis (H3 digit 1) has no + * child cell (the deleted subsequence). To handle pentagonal input, we build + * the digit sequence as if it were a hexagon, then skip over the deleted + * subsequences (edges that land on the missing K-axis side). + */ + +#include +#include + +#include "coordijk.h" +#include "h3Assert.h" +#include "h3Index.h" +#include "iterators.h" +#include "mathExtensions.h" + +// H3 digit at each of the WALK_CYCLE_LEN boundary walk positions. +// The 6 non-center digits {1,5,4,6,2,3} in counter-clockwise order around +// the hexagon, each repeated 3x for the fractal subdivision at each side. +// Note: walk_digit[i] == edge_dir[(i/3 + 1) % 6]; the same cyclic sequence +// of digits, rotated by 1 and each element repeated 3x. +static const Direction walk_digit[] = { + K_AXES_DIGIT, K_AXES_DIGIT, K_AXES_DIGIT, // 1,1,1 + IK_AXES_DIGIT, IK_AXES_DIGIT, IK_AXES_DIGIT, // 5,5,5 + I_AXES_DIGIT, I_AXES_DIGIT, I_AXES_DIGIT, // 4,4,4 + IJ_AXES_DIGIT, IJ_AXES_DIGIT, IJ_AXES_DIGIT, // 6,6,6 + J_AXES_DIGIT, J_AXES_DIGIT, J_AXES_DIGIT, // 2,2,2 + JK_AXES_DIGIT, JK_AXES_DIGIT, JK_AXES_DIGIT}; // 3,3,3 + +#define WALK_CYCLE_LEN ((int)(sizeof(walk_digit) / sizeof(walk_digit[0]))) + +// H3 edge direction at each edge index, in counter-clockwise order around the +// hexagon. Stored in the directed edge's reserved bits. +// Numeric values: {3, 1, 5, 4, 6, 2} +static const Direction edge_dir[] = {JK_AXES_DIGIT, K_AXES_DIGIT, + IK_AXES_DIGIT, I_AXES_DIGIT, + IJ_AXES_DIGIT, J_AXES_DIGIT}; + +/** + * Advance the walk along origin cells on the Gosper island boundary, + * updating the child digits. + * + * Each resolution level cycles through WALK_CYCLE_LEN walk positions (6 sides * + * 3 segments per side). Across resolutions, each step at a coarser level maps + * to 3 steps at the next finer level. When the coarser level moves to the next + * cell, the finer level shifts back by 6 positions. + * + * The point where each level moves to the next cell depends on + * resolution parity: Class II (even r) at walkPos % 3 == 0, + * Class III (odd r) at walkPos % 3 == 1. + * + * @return true if the origin cell changed. + */ +static bool advanceOriginCell(int8_t *walkPos, H3Index *h, int8_t r, + int8_t parentRes) { + int prevDigit = H3_GET_INDEX_DIGIT(*h, r); + + // When moving to the next cell, recurse to advance the parent resolution. + if ((r > parentRes + 1) && (walkPos[r] % 3 == r % 2)) { + bool parentChanged = advanceOriginCell(walkPos, h, r - 1, parentRes); + if (parentChanged) { + // New parent cell: shift back by 6 positions. + walkPos[r] -= 6; + } + } + + // Advance to next position: Increase the cyclic index by one, but + // we use + (WALK_CYCLE_LEN + 1) == +1 (mod WALK_CYCLE_LEN), so that it + // stays positive, even after potentially subtracting by 6 above. + walkPos[r] = (walkPos[r] + WALK_CYCLE_LEN + 1) % WALK_CYCLE_LEN; + + // Update the child digit + Direction newDigit = walk_digit[walkPos[r]]; + H3_SET_INDEX_DIGIT(*h, r, newDigit); + + // A change in the finest digit is sufficient to detect a cell change. + return newDigit != prevDigit; +} + +/** + * Advance to the next boundary edge (origin cell + edge direction). + */ +static void advanceEdge(IterEdgesGosper *iter) { + bool cellChanged = advanceOriginCell(iter->_walkPos, &(iter->e), + iter->_childRes, iter->_parentRes); + + // If we move to a new origin cell, shift the _edgePos back by 2 + if (cellChanged) { + iter->_edgePos -= 2; + } + + // Advance to next edge. Increase the cyclic index by one, but + // we use +7 == +1 (mod 6) so that it stays positive, + // even after potentially subtracting by 2 above. + iter->_edgePos = (iter->_edgePos + 7) % 6; + H3_SET_RESERVED_BITS(iter->e, edge_dir[iter->_edgePos]); +} + +// Skip edges corresponding to deleted subsequences of a pentagon. +static void skipPentagonEdges(IterEdgesGosper *iter) { + while (iter->_isPentagon && ALWAYS(iter->e != H3_NULL) && + H3_GET_INDEX_DIGIT(iter->e, iter->_parentRes + 1) == 1) { + advanceEdge(iter); + } +} + +void iterStepGosper(IterEdgesGosper *iter) { + if (iter->e == H3_NULL) return; + + iter->remaining--; + if (iter->remaining <= 0) { + iter->e = H3_NULL; + return; + } + + if (iter->_parentRes == iter->_childRes) { + // Easy case of 0 resolution delta: + // Cycle through edge indices, skipping the + // deleted subsequence for pentagons. + iter->_edgePos++; + if (iter->_isPentagon && edge_dir[iter->_edgePos] == K_AXES_DIGIT) { + iter->_edgePos++; + } + + H3_SET_RESERVED_BITS(iter->e, edge_dir[iter->_edgePos]); + } else { + // Resolution delta > 0 + advanceEdge(iter); + + // Skip deleted subsequences for pentagons. + skipPentagonEdges(iter); + } +} + +// Note: No input validation here. This is an internal function where +// validation is expected to be done upstream. If this is ever added to the +// public API, we'll need to validate h, childRes >= parentRes, +// and childRes <= MAX_H3_RES. +IterEdgesGosper iterInitGosper(H3Index h, int childRes) { + int parentRes = H3_GET_RESOLUTION(h); + bool isPent = H3_EXPORT(isPentagon)(h); + + IterEdgesGosper iter = { + .e = h, + .remaining = 0, + + ._edgePos = 0, + ._walkPos = {0}, + + ._parentRes = parentRes, + ._childRes = childRes, + ._isPentagon = isPent, + }; + + // Set number of edges for the iterator to yield + int numSides = isPent ? 5 : 6; + iter.remaining = numSides * _ipow(3, childRes - parentRes); + + // Same-resolution fast path: no digit/resolution setup needed + if (parentRes == childRes) { + H3_SET_MODE(iter.e, H3_DIRECTEDEDGE_MODE); + // Skip the deleted subsequence for pentagons (if edge_dir happens + // to start there) + if (NEVER(isPent && edge_dir[iter._edgePos] == K_AXES_DIGIT)) { + iter._edgePos++; + } + H3_SET_RESERVED_BITS(iter.e, edge_dir[iter._edgePos]); + return iter; + } + + // Set resolution to child level and initialize boundary walk. + // The first child level starts at walk position 0. + // Subsequent levels start at either 14 or 16, depending on + // if those resolutions are Class II or III. + H3_SET_RESOLUTION(iter.e, childRes); + for (int r = parentRes + 1; r <= childRes; r++) { + if (r == parentRes + 1) { + iter._walkPos[r] = 0; + } else { + iter._walkPos[r] = (r % 2) ? 14 : 16; + } + H3_SET_INDEX_DIGIT(iter.e, r, walk_digit[iter._walkPos[r]]); + } + + // Set edge mode and initial direction + H3_SET_MODE(iter.e, H3_DIRECTEDEDGE_MODE); + H3_SET_RESERVED_BITS(iter.e, edge_dir[iter._edgePos]); + + // For pentagons, the walk may start on a deleted subsequence. + // Advance past those positions to reach the first valid edge. + skipPentagonEdges(&iter); + + return iter; +}