From df1b842fcd484315d4ca441a36634d82a42b14b3 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 26 Nov 2025 23:46:34 -0500 Subject: [PATCH 01/74] area.c --- justfile | 45 +++++++++++++++++++++++++++ src/h3lib/lib/area.c | 74 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 119 insertions(+) create mode 100644 justfile create mode 100644 src/h3lib/lib/area.c diff --git a/justfile b/justfile new file mode 100644 index 0000000000..8ab59e9ced --- /dev/null +++ b/justfile @@ -0,0 +1,45 @@ +# TODO: Remove before landing + +init: purge + mkdir build + +build: + cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + +profile: init + cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make + # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + xcrun xctrace record \ + --template 'Time Profiler' \ + --output h3-prof.trace \ + --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon + open h3-prof.trace + +purge: + rm -rf build + rm -rf *.trace + rm -rf .ipynb_checkpoints + +test: build + # ./build/bin/testH3CellAreaExhaustive + # ./build/bin/testEdgeCellsToPoly + ./build/bin/testDirectedEdge + # ./build/bin/testArea + # just test-slow + +time: + time ./build/bin/testArea + +test-fast: build + cd build; make test-fast + +test-slow: build + cd build; make test + +bench: build + ./build/bin/benchmarkCellsToLinkedMultiPolygon + # ./build/bin/benchmarkCellsToMultiPolygon + # ./build/bin/benchmarkDirectedEdge + +fuzz: + cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c new file mode 100644 index 0000000000..e14677f307 --- /dev/null +++ b/src/h3lib/lib/area.c @@ -0,0 +1,74 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "adder.h" +#include "alloc.h" +#include "h3api.h" + + +void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop) { + H3_MEMORY(free)(loop->verts); + loop->verts = NULL; +} + +static double cagnoli(LatLng x, LatLng y) { + x.lat = x.lat / 2.0 + M_PI / 4.0; + y.lat = y.lat / 2.0 + M_PI / 4.0; + + double sa = sin(x.lat) * sin(y.lat); + double ca = cos(x.lat) * cos(y.lat); + + double d = y.lng - x.lng; + double sd = sin(d); + double cd = cos(d); + + return -2.0 * atan2(sa * sd, sa * cd + ca); +} + +static double vertArea(LatLng *verts, int numVerts) { + Kahan k = {0.0, 0.0}; + + for (int i = 0; i < numVerts; i++) { + int j = (i + 1) % numVerts; + kadd(&k, cagnoli(verts[i], verts[j])); + } + + if (k.sum < 0.0) { + kadd(&k, 4.0 * M_PI); + } + return k.sum; +} + +H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { + *out = vertArea(loop.verts, loop.numVerts); + return E_SUCCESS; +} + +/** + * Area of H3 cell in radians^2. + * + * The area is calculated by breaking the cell into spherical triangles and + * summing up their areas. Note that some H3 cells (hexagons and pentagons) + * are irregular, and have more than 6 or 5 sides. + * + * todo: optimize the computation by re-using the edges shared between triangles + * + * @param cell H3 cell + * @param out cell area in radians^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { + CellBoundary cb; + H3Error err = H3_EXPORT(cellToBoundary)(cell, &cb); + if (err) { + return err; + } + *out = vertArea(cb.verts, cb.numVerts); + + return E_SUCCESS; +} From b6d9ddb9c16ba6881f8d76f9c79329fb881407da Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 26 Nov 2025 23:58:13 -0500 Subject: [PATCH 02/74] adding in more stuff --- CMakeLists.txt | 5 ++- src/apps/testapps/testArea.c | 56 +++++++++++++++++++++++++++ src/h3lib/include/adder.h | 16 ++++++++ src/h3lib/include/h3api.h.in | 3 ++ src/h3lib/lib/latLng.c | 74 ------------------------------------ 5 files changed, 79 insertions(+), 75 deletions(-) create mode 100644 src/apps/testapps/testArea.c create mode 100644 src/h3lib/include/adder.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 7eb6dbddb5..828c02bcef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -170,6 +170,7 @@ set(LIB_SOURCE_FILES src/h3lib/include/constants.h src/h3lib/include/coordijk.h src/h3lib/include/algos.h + src/h3lib/include/adder.h src/h3lib/lib/h3Assert.c src/h3lib/lib/algos.c src/h3lib/lib/coordijk.c @@ -188,7 +189,8 @@ set(LIB_SOURCE_FILES src/h3lib/lib/iterators.c src/h3lib/lib/vertexGraph.c src/h3lib/lib/faceijk.c - src/h3lib/lib/baseCells.c) + src/h3lib/lib/baseCells.c + src/h3lib/lib/area.c) set(APP_SOURCE_FILES src/apps/applib/include/kml.h src/apps/applib/include/benchmark.h @@ -278,6 +280,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c new file mode 100644 index 0000000000..1482e74176 --- /dev/null +++ b/src/apps/testapps/testArea.c @@ -0,0 +1,56 @@ +#include +#include +#include +#include + +#include "h3api.h" +#include "iterators.h" +#include "test.h" +#include "utility.h" +#include "adder.h" + + +static void do_res_sum(int res) { + Kahan k = {0,0}; + + for (IterCellsResolution iter = iterInitRes(res); iter.h; + iterStepRes(&iter)) { + double cellArea; + H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); + + kadd(&k, cellArea); + } + + double area = k.sum; + double diff = fabs(area - 4 * M_PI); + printf("res: %d, diff: %e\n", res, diff); +} + +SUITE(test_for_area) { + TEST(some_area_test) { + // the numerical test is how close to 4*pi do we get when adding up all + // cells at smaller and smaller resolution? + + printf("\n"); + + do_res_sum(0); + do_res_sum(1); + do_res_sum(2); + do_res_sum(3); + do_res_sum(4); + // do_res_sum(5); + // do_res_sum(6); + // do_res_sum(7); + // do_res_sum(8); + + // res: 0, diff: 7.105427e-15 + // res: 1, diff: 3.907985e-14 + // res: 2, diff: 1.421085e-13 + // res: 3, diff: 4.689582e-13 + // res: 4, diff: 2.161826e-12 + // res: 5, diff: 1.543654e-12 + // res: 6, diff: 1.768363e-11 + // res: 7, diff: 1.490719e-11 + // res: 8, diff: 2.937917e-11 + } +} diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h new file mode 100644 index 0000000000..60e19891d6 --- /dev/null +++ b/src/h3lib/include/adder.h @@ -0,0 +1,16 @@ +#ifndef ADDER_H +#define ADDER_H + +typedef struct { + double sum; // running total + double c; // compensation term +} Kahan; + +static inline void kadd(Kahan *k, double x) { + double y = x - k->c; + double t = k->sum + y; + k->c = (t - k->sum) - y; + k->sum = t; +} + +#endif diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index dfca733620..a3159cffee 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,6 +348,9 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ +DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop); +DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) + /** @defgroup degsToRads degsToRads * Functions for degsToRads * @{ diff --git a/src/h3lib/lib/latLng.c b/src/h3lib/lib/latLng.c index f90d9ee2b4..13bf196e8e 100644 --- a/src/h3lib/lib/latLng.c +++ b/src/h3lib/lib/latLng.c @@ -353,80 +353,6 @@ H3Error H3_EXPORT(getNumCells)(int res, int64_t *out) { return E_SUCCESS; } -/** - * Surface area in radians^2 of spherical triangle on unit sphere. - * - * For the math, see: - * https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess - * - * @param a length of triangle side A in radians - * @param b length of triangle side B in radians - * @param c length of triangle side C in radians - * - * @return area in radians^2 of triangle on unit sphere - */ -double triangleEdgeLengthsToArea(double a, double b, double c) { - double s = (a + b + c) * 0.5; - - a = (s - a) * 0.5; - b = (s - b) * 0.5; - c = (s - c) * 0.5; - s = s * 0.5; - - return 4 * atan(sqrt(tan(s) * tan(a) * tan(b) * tan(c))); -} - -/** - * Compute area in radians^2 of a spherical triangle, given its vertices. - * - * @param a vertex lat/lng in radians - * @param b vertex lat/lng in radians - * @param c vertex lat/lng in radians - * - * @return area of triangle on unit sphere, in radians^2 - */ -double triangleArea(const LatLng *a, const LatLng *b, const LatLng *c) { - return triangleEdgeLengthsToArea(H3_EXPORT(greatCircleDistanceRads)(a, b), - H3_EXPORT(greatCircleDistanceRads)(b, c), - H3_EXPORT(greatCircleDistanceRads)(c, a)); -} - -/** - * Area of H3 cell in radians^2. - * - * The area is calculated by breaking the cell into spherical triangles and - * summing up their areas. Note that some H3 cells (hexagons and pentagons) - * are irregular, and have more than 6 or 5 sides. - * - * todo: optimize the computation by re-using the edges shared between triangles - * - * @param cell H3 cell - * @param out cell area in radians^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { - LatLng c; - CellBoundary cb; - H3Error err = H3_EXPORT(cellToLatLng)(cell, &c); - if (err) { - return err; - } - err = H3_EXPORT(cellToBoundary)(cell, &cb); - if (NEVER(err)) { - // Uncoverable because cellToLatLng will have returned an error already - return err; - } - - double area = 0.0; - for (int i = 0; i < cb.numVerts; i++) { - int j = (i + 1) % cb.numVerts; - area += triangleArea(&cb.verts[i], &cb.verts[j], &c); - } - - *out = area; - return E_SUCCESS; -} - /** * Area of H3 cell in kilometers^2. * From 8080e568e290766f61c0f049fa762c332f08a705 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 00:03:06 -0500 Subject: [PATCH 03/74] merp --- CMakeTests.cmake | 1 + justfile | 4 ++-- src/apps/testapps/testArea.c | 5 ++--- src/h3lib/include/h3api.h.in | 4 ++-- src/h3lib/lib/area.c | 3 +-- 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 1901a61fac..cb705fc94a 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,6 +259,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(testArea src/apps/testapps/testArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index 8ab59e9ced..c9992e3ce8 100644 --- a/justfile +++ b/justfile @@ -23,8 +23,8 @@ purge: test: build # ./build/bin/testH3CellAreaExhaustive # ./build/bin/testEdgeCellsToPoly - ./build/bin/testDirectedEdge - # ./build/bin/testArea + # ./build/bin/testDirectedEdge + ./build/bin/testArea # just test-slow time: diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c index 1482e74176..b103c950ec 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/testapps/testArea.c @@ -3,15 +3,14 @@ #include #include +#include "adder.h" #include "h3api.h" #include "iterators.h" #include "test.h" #include "utility.h" -#include "adder.h" - static void do_res_sum(int res) { - Kahan k = {0,0}; + Kahan k = {0, 0}; for (IterCellsResolution iter = iterInitRes(res); iter.h; iterStepRes(&iter)) { diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index a3159cffee..042b778471 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,8 +348,8 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop); -DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) +DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop); +DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); /** @defgroup degsToRads degsToRads * Functions for degsToRads diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index e14677f307..08a227629c 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -10,8 +10,7 @@ #include "alloc.h" #include "h3api.h" - -void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop) { +void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { H3_MEMORY(free)(loop->verts); loop->verts = NULL; } From 5c8950ade3992c4c55a56e3f00261645917c1c33 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 00:06:05 -0500 Subject: [PATCH 04/74] kahan accumulator --- src/apps/testapps/testH3CellAreaExhaustive.c | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 13261dae15..c0bb793de7 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -136,15 +136,16 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - double area = 0.0; - for (IterCellsResolution iter = iterInitRes(res); iter.h; - iterStepRes(&iter)) { + Kahan k = {0.0, 0.0}; + IterCellsResolution iter = iterInitRes(res); + + for (; iter.h; iterStepRes(&iter)) { double cellArea; t_assertSuccess((*cell_area)(iter.h, &cellArea)); - area += cellArea; + kadd(&k, cellArea); } - t_assert(fabs(area - target) < tol, + t_assert(fabs(k.sum - target) < tol, "sum of all cells should give earth area"); } From b82cab22a46877986dd70c6d9ab58e967ac54ed8 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 11:18:35 -0500 Subject: [PATCH 05/74] neumaier_add --- src/apps/testapps/testArea.c | 15 ++++++--- src/apps/testapps/testH3CellAreaExhaustive.c | 1 + src/h3lib/include/adder.h | 32 +++++++++++++++++++- src/h3lib/lib/area.c | 2 +- 4 files changed, 44 insertions(+), 6 deletions(-) diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c index b103c950ec..813b542d14 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/testapps/testArea.c @@ -9,6 +9,13 @@ #include "test.h" #include "utility.h" +// TODO: turn this into a benchmark, report accuracy and time +// see if we get any benefit by avoiding the repeated sin/cos computation + +// TODO: document algos: https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements +// TODO: better more general name for accumulator struct and functions +// TOOD: demonstrate the area alg works for global polygons. + static void do_res_sum(int res) { Kahan k = {0, 0}; @@ -20,7 +27,7 @@ static void do_res_sum(int res) { kadd(&k, cellArea); } - double area = k.sum; + double area = kresult(k); double diff = fabs(area - 4 * M_PI); printf("res: %d, diff: %e\n", res, diff); } @@ -37,9 +44,9 @@ SUITE(test_for_area) { do_res_sum(2); do_res_sum(3); do_res_sum(4); - // do_res_sum(5); - // do_res_sum(6); - // do_res_sum(7); + do_res_sum(5); + do_res_sum(6); + do_res_sum(7); // do_res_sum(8); // res: 0, diff: 7.105427e-15 diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index c0bb793de7..3f39fac7bb 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -23,6 +23,7 @@ #include +#include "adder.h" #include "constants.h" #include "h3api.h" #include "iterators.h" diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 60e19891d6..a9763be3ed 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -6,11 +6,41 @@ typedef struct { double c; // compensation term } Kahan; -static inline void kadd(Kahan *k, double x) { +static inline void simple_add(Kahan *k, double x) { k->sum += x; } + +static inline void kahan_add(Kahan *k, double x) { double y = x - k->c; double t = k->sum + y; k->c = (t - k->sum) - y; k->sum = t; } +static inline void neumaier_add(Kahan *k, double x) { + double t = k->sum + x; + + if (fabs(k->sum) >= fabs(x)) { + k->c += (k->sum - t) + x; + } else { + k->c += (x - t) + k->sum; + } + + k->sum = t; +} + +static inline double simple_result(Kahan k) { return k.sum; } +static inline double kahan_result(Kahan k) { return k.sum; } +static inline double neumaier_result(Kahan k) { return k.sum + k.c; } + +static inline void kadd(Kahan *k, double x) { + // simple_add(k, x); + // kahan_add(k, x); + neumaier_add(k, x); +} + +static inline double kresult(Kahan k) { + // return simple_result(k); + // return kahan_result(k); + return neumaier_result(k); +} + #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 08a227629c..1eef4940e9 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -40,7 +40,7 @@ static double vertArea(LatLng *verts, int numVerts) { if (k.sum < 0.0) { kadd(&k, 4.0 * M_PI); } - return k.sum; + return kresult(k); } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { From ec26cb85f39618d649203399721e72fc5a8d0882 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 11:41:43 -0500 Subject: [PATCH 06/74] notes --- src/apps/testapps/testArea.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c index 813b542d14..d927210ff8 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/testapps/testArea.c @@ -12,15 +12,18 @@ // TODO: turn this into a benchmark, report accuracy and time // see if we get any benefit by avoiding the repeated sin/cos computation -// TODO: document algos: https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements +// TODO: document algos: +// https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements // TODO: better more general name for accumulator struct and functions // TOOD: demonstrate the area alg works for global polygons. +// TODO: is kadd faster without pointers? k = kadd(k, cellArea); + static void do_res_sum(int res) { Kahan k = {0, 0}; + IterCellsResolution iter = iterInitRes(res); - for (IterCellsResolution iter = iterInitRes(res); iter.h; - iterStepRes(&iter)) { + for (; iter.h; iterStepRes(&iter)) { double cellArea; H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); From 74e277ff43312452ab8b6d0eaa922a2d38f4b579 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 11:52:27 -0500 Subject: [PATCH 07/74] benchmark area --- CMakeLists.txt | 6 ++-- CMakeTests.cmake | 1 - justfile | 6 ++-- .../testArea.c => benchmarks/benchmarkArea.c} | 35 +++++-------------- 4 files changed, 16 insertions(+), 32 deletions(-) rename src/apps/{testapps/testArea.c => benchmarks/benchmarkArea.c} (56%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 828c02bcef..7e4f88f7d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -280,7 +280,6 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c - src/apps/testapps/testArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c @@ -320,7 +319,8 @@ set(OTHER_SOURCE_FILES src/apps/benchmarks/benchmarkDirectedEdge.c src/apps/benchmarks/benchmarkVertex.c src/apps/benchmarks/benchmarkIsValidCell.c - src/apps/benchmarks/benchmarkH3Api.c) + src/apps/benchmarks/benchmarkH3Api.c + src/apps/benchmarks/benchmarkArea.c) set(ALL_SOURCE_FILES ${LIB_SOURCE_FILES} ${APP_SOURCE_FILES} ${TEST_APP_SOURCE_FILES} @@ -675,6 +675,8 @@ if(BUILD_BENCHMARKS) src/apps/benchmarks/benchmarkPolygonToCells.c) add_h3_benchmark(benchmarkPolygonToCellsExperimental src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c) + add_h3_benchmark(benchmarkArea + src/apps/benchmarks/benchmarkArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_benchmark(benchmarkPolygon src/apps/benchmarks/benchmarkPolygon.c) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index cb705fc94a..1901a61fac 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,7 +259,6 @@ 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(testArea src/apps/testapps/testArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index c9992e3ce8..5d7ff258c8 100644 --- a/justfile +++ b/justfile @@ -24,8 +24,9 @@ test: build # ./build/bin/testH3CellAreaExhaustive # ./build/bin/testEdgeCellsToPoly # ./build/bin/testDirectedEdge - ./build/bin/testArea + # ./build/bin/testArea # just test-slow + just test-fast time: time ./build/bin/testArea @@ -37,7 +38,8 @@ test-slow: build cd build; make test bench: build - ./build/bin/benchmarkCellsToLinkedMultiPolygon + ./build/bin/benchmarkArea + # ./build/bin/benchmarkCellsToLinkedMultiPolygon # ./build/bin/benchmarkCellsToMultiPolygon # ./build/bin/benchmarkDirectedEdge diff --git a/src/apps/testapps/testArea.c b/src/apps/benchmarks/benchmarkArea.c similarity index 56% rename from src/apps/testapps/testArea.c rename to src/apps/benchmarks/benchmarkArea.c index d927210ff8..1ffe6fa858 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -4,9 +4,9 @@ #include #include "adder.h" +#include "benchmark.h" #include "h3api.h" #include "iterators.h" -#include "test.h" #include "utility.h" // TODO: turn this into a benchmark, report accuracy and time @@ -35,31 +35,12 @@ static void do_res_sum(int res) { printf("res: %d, diff: %e\n", res, diff); } -SUITE(test_for_area) { - TEST(some_area_test) { - // the numerical test is how close to 4*pi do we get when adding up all - // cells at smaller and smaller resolution? +BEGIN_BENCHMARKS(); - printf("\n"); - - do_res_sum(0); - do_res_sum(1); - do_res_sum(2); - do_res_sum(3); - do_res_sum(4); - do_res_sum(5); - do_res_sum(6); - do_res_sum(7); - // do_res_sum(8); - - // res: 0, diff: 7.105427e-15 - // res: 1, diff: 3.907985e-14 - // res: 2, diff: 1.421085e-13 - // res: 3, diff: 4.689582e-13 - // res: 4, diff: 2.161826e-12 - // res: 5, diff: 1.543654e-12 - // res: 6, diff: 1.768363e-11 - // res: 7, diff: 1.490719e-11 - // res: 8, diff: 2.937917e-11 +BENCHMARK(directedEdgeToBoundary, 1, { + for (int i = 0; i < 7; i++) { + do_res_sum(i); } -} +}); + +END_BENCHMARKS(); From c915b5183011839bb211a305241454134ee8c759 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 22:11:55 -0500 Subject: [PATCH 08/74] turn it to 11. --- src/apps/benchmarks/benchmarkArea.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 1ffe6fa858..ae006a09ef 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -38,7 +38,7 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - for (int i = 0; i < 7; i++) { + for (int i = 0; i <= 8; i++) { do_res_sum(i); } }); From a74b9d949042e56873150cf71d3a9af828544386 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 23:31:41 -0500 Subject: [PATCH 09/74] sadly, no improvement from reusing some trig --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/h3lib/lib/area.c | 54 ++++++++++++++++++++++++++--- 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index ae006a09ef..d246493941 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -38,7 +38,7 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - for (int i = 0; i <= 8; i++) { + for (int i = 0; i <= 6; i++) { do_res_sum(i); } }); diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 1eef4940e9..bd8b5551b5 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -15,7 +15,53 @@ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { loop->verts = NULL; } -static double cagnoli(LatLng x, LatLng y) { +typedef struct { + double s; + double c; + double lng; +} LatLngPre; + +static inline LatLngPre precompute_latlng(LatLng x) { + double lat = x.lat / 2.0 + M_PI / 4.0; + + LatLngPre llp = {.s = sin(lat), .c = cos(lat), .lng = x.lng}; + + return llp; +} + +static inline double cagnoli_pre(LatLngPre x, LatLngPre y) { + double sa = x.s * y.s; + double ca = x.c * y.c; + + double d = y.lng - x.lng; + double sd = sin(d); + double cd = cos(d); + + return -2.0 * atan2(sa * sd, sa * cd + ca); +} + +static inline double vertArea2(LatLng *verts, int numVerts) { + Kahan k = {0.0, 0.0}; + + LatLngPre x = precompute_latlng(verts[0]); + LatLngPre first = x; + LatLngPre y; + + for (int i = 1; i < numVerts; i++) { + y = precompute_latlng(verts[i]); + + kadd(&k, cagnoli_pre(x, y)); + x = y; + } + kadd(&k, cagnoli_pre(x, first)); + + if (k.sum < 0.0) { + kadd(&k, 4.0 * M_PI); + } + return kresult(k); +} + +static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -29,7 +75,7 @@ static double cagnoli(LatLng x, LatLng y) { return -2.0 * atan2(sa * sd, sa * cd + ca); } -static double vertArea(LatLng *verts, int numVerts) { +static inline double vertArea(LatLng *verts, int numVerts) { Kahan k = {0.0, 0.0}; for (int i = 0; i < numVerts; i++) { @@ -44,7 +90,7 @@ static double vertArea(LatLng *verts, int numVerts) { } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - *out = vertArea(loop.verts, loop.numVerts); + *out = vertArea2(loop.verts, loop.numVerts); return E_SUCCESS; } @@ -67,7 +113,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { if (err) { return err; } - *out = vertArea(cb.verts, cb.numVerts); + *out = vertArea2(cb.verts, cb.numVerts); return E_SUCCESS; } From 98ecbee87cb17e308beb549854db55bdf1c5744e Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 23:32:43 -0500 Subject: [PATCH 10/74] back to original --- src/h3lib/lib/area.c | 50 ++------------------------------------------ 1 file changed, 2 insertions(+), 48 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index bd8b5551b5..ec7c39ecc4 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -15,52 +15,6 @@ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { loop->verts = NULL; } -typedef struct { - double s; - double c; - double lng; -} LatLngPre; - -static inline LatLngPre precompute_latlng(LatLng x) { - double lat = x.lat / 2.0 + M_PI / 4.0; - - LatLngPre llp = {.s = sin(lat), .c = cos(lat), .lng = x.lng}; - - return llp; -} - -static inline double cagnoli_pre(LatLngPre x, LatLngPre y) { - double sa = x.s * y.s; - double ca = x.c * y.c; - - double d = y.lng - x.lng; - double sd = sin(d); - double cd = cos(d); - - return -2.0 * atan2(sa * sd, sa * cd + ca); -} - -static inline double vertArea2(LatLng *verts, int numVerts) { - Kahan k = {0.0, 0.0}; - - LatLngPre x = precompute_latlng(verts[0]); - LatLngPre first = x; - LatLngPre y; - - for (int i = 1; i < numVerts; i++) { - y = precompute_latlng(verts[i]); - - kadd(&k, cagnoli_pre(x, y)); - x = y; - } - kadd(&k, cagnoli_pre(x, first)); - - if (k.sum < 0.0) { - kadd(&k, 4.0 * M_PI); - } - return kresult(k); -} - static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -90,7 +44,7 @@ static inline double vertArea(LatLng *verts, int numVerts) { } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - *out = vertArea2(loop.verts, loop.numVerts); + *out = vertArea(loop.verts, loop.numVerts); return E_SUCCESS; } @@ -113,7 +67,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { if (err) { return err; } - *out = vertArea2(cb.verts, cb.numVerts); + *out = vertArea(cb.verts, cb.numVerts); return E_SUCCESS; } From db21eb0bc0ba056a474f873c038dc8145546a5ff Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 00:12:41 -0500 Subject: [PATCH 11/74] simplify --- src/h3lib/lib/area.c | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index ec7c39ecc4..decc096afd 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -10,11 +10,6 @@ #include "alloc.h" #include "h3api.h" -void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { - H3_MEMORY(free)(loop->verts); - loop->verts = NULL; -} - static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -29,22 +24,19 @@ static inline double cagnoli(LatLng x, LatLng y) { return -2.0 * atan2(sa * sd, sa * cd + ca); } -static inline double vertArea(LatLng *verts, int numVerts) { +H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { Kahan k = {0.0, 0.0}; - for (int i = 0; i < numVerts; i++) { - int j = (i + 1) % numVerts; - kadd(&k, cagnoli(verts[i], verts[j])); + for (int i = 0; i < loop.numVerts; i++) { + int j = (i + 1) % loop.numVerts; + kadd(&k, cagnoli(loop.verts[i], loop.verts[j])); } if (k.sum < 0.0) { kadd(&k, 4.0 * M_PI); } - return kresult(k); -} -H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - *out = vertArea(loop.verts, loop.numVerts); + *out = kresult(k); return E_SUCCESS; } @@ -67,7 +59,17 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { if (err) { return err; } - *out = vertArea(cb.verts, cb.numVerts); + + GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; + err = H3_EXPORT(geoLoopArea)(loop, out); + if (err) { + return err; + } return E_SUCCESS; } + +void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { + H3_MEMORY(free)(loop->verts); + loop->verts = NULL; +} From 04640ef8b648f90e8cdfcafce3ee7e80e770fbb5 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 20:49:25 -0500 Subject: [PATCH 12/74] try adding constants --- src/h3lib/lib/area.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index decc096afd..cec3c90ddd 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -9,6 +9,7 @@ #include "adder.h" #include "alloc.h" #include "h3api.h" +#include "constants.h" static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; From 2349ab4c333a35b9bd8aff560be60b0a3933e949 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 20:53:27 -0500 Subject: [PATCH 13/74] format, we must --- src/h3lib/lib/area.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index cec3c90ddd..203f039c9a 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -8,8 +8,8 @@ #include "adder.h" #include "alloc.h" -#include "h3api.h" #include "constants.h" +#include "h3api.h" static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; From 0dd6a5eec25bfb3a3ff93f934a9d52c36c4b11a7 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 22:10:07 -0500 Subject: [PATCH 14/74] Adder docs --- src/apps/benchmarks/benchmarkArea.c | 17 ++--- src/apps/testapps/testH3CellAreaExhaustive.c | 6 +- src/h3lib/include/adder.h | 78 ++++++++++++++------ src/h3lib/lib/area.c | 10 +-- 4 files changed, 67 insertions(+), 44 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index d246493941..69ce545cfb 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,29 +9,22 @@ #include "iterators.h" #include "utility.h" -// TODO: turn this into a benchmark, report accuracy and time -// see if we get any benefit by avoiding the repeated sin/cos computation - -// TODO: document algos: -// https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements -// TODO: better more general name for accumulator struct and functions // TOOD: demonstrate the area alg works for global polygons. // TODO: is kadd faster without pointers? k = kadd(k, cellArea); static void do_res_sum(int res) { - Kahan k = {0, 0}; + Adder adder = {0, 0}; + double cellArea; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { - double cellArea; H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); - - kadd(&k, cellArea); + kadd(&adder, cellArea); } - double area = kresult(k); - double diff = fabs(area - 4 * M_PI); + double total_area = kresult(adder); + double diff = fabs(total_area - 4 * M_PI); printf("res: %d, diff: %e\n", res, diff); } diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 3f39fac7bb..b464014b46 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -137,16 +137,16 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - Kahan k = {0.0, 0.0}; + Adder adder = {0.0, 0.0}; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { double cellArea; t_assertSuccess((*cell_area)(iter.h, &cellArea)); - kadd(&k, cellArea); + kadd(&adder, cellArea); } - t_assert(fabs(k.sum - target) < tol, + t_assert(fabs(kresult(adder) - target) < tol, "sum of all cells should give earth area"); } diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index a9763be3ed..e391987964 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -1,46 +1,76 @@ #ifndef ADDER_H #define ADDER_H +/* +This is a header-only implementation of a "compensated summation" algorithm, +which allows us to add up sequences of floating point numbers with better +precision than naive summation, especially when the terms in the sum +can vary in magnitude. +See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm + +This is useful when computing the area of (multi-)polygons, which can involve +adding many small terms to large aggregates. +For example, D3 uses this when computing polygonal areas via the `Kahan` class: +https://github.com/d3/d3-geo/blob/main/src/area.js + +There are a few potential algorithms we might consider for summation: + +1. Naive sum +2. Kahan summation +3. Neumaier summation +4. Other approaches like pairwise summation, or Pythons `fsum` + +We considered the first three for simplicity, and settled on Kahan summation, +as it acheives noticably better accuracy than naive summation and almost as good +accuracy as Neumaier, while only being slightly slower than naive, but slightly +faster and simpler than Neumaier. + +See also: https://github.com/python/cpython/issues/100425 for discussions of +tradeoffs between Kahan, Neumaier, and `fsum`. +*/ + typedef struct { double sum; // running total double c; // compensation term -} Kahan; +} Adder; -static inline void simple_add(Kahan *k, double x) { k->sum += x; } +static inline void simple_add(Adder *adder, double x) { adder->sum += x; } -static inline void kahan_add(Kahan *k, double x) { - double y = x - k->c; - double t = k->sum + y; - k->c = (t - k->sum) - y; - k->sum = t; +static inline void kahan_add(Adder *adder, double x) { + double y = x - adder->c; + double t = adder->sum + y; + adder->c = (t - adder->sum) - y; + adder->sum = t; } -static inline void neumaier_add(Kahan *k, double x) { - double t = k->sum + x; +static inline void neumaier_add(Adder *adder, double x) { + double t = adder->sum + x; - if (fabs(k->sum) >= fabs(x)) { - k->c += (k->sum - t) + x; + if (fabs(adder->sum) >= fabs(x)) { + adder->c += (adder->sum - t) + x; } else { - k->c += (x - t) + k->sum; + adder->c += (x - t) + adder->sum; } - k->sum = t; + adder->sum = t; } -static inline double simple_result(Kahan k) { return k.sum; } -static inline double kahan_result(Kahan k) { return k.sum; } -static inline double neumaier_result(Kahan k) { return k.sum + k.c; } +static inline double simple_result(Adder adder) { return adder.sum; } +static inline double kahan_result(Adder adder) { return adder.sum; } +static inline double neumaier_result(Adder adder) { + return adder.sum + adder.c; +} -static inline void kadd(Kahan *k, double x) { - // simple_add(k, x); - // kahan_add(k, x); - neumaier_add(k, x); +static inline void kadd(Adder *adder, double x) { + // simple_add(adder, x); + // kahan_add(adder, x); + neumaier_add(adder, x); } -static inline double kresult(Kahan k) { - // return simple_result(k); - // return kahan_result(k); - return neumaier_result(k); +static inline double kresult(Adder adder) { + // return simple_result(adder); + // return kahan_result(adder); + return neumaier_result(adder); } #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 203f039c9a..a10e337dcf 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -26,18 +26,18 @@ static inline double cagnoli(LatLng x, LatLng y) { } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - Kahan k = {0.0, 0.0}; + Adder adder = {0.0, 0.0}; for (int i = 0; i < loop.numVerts; i++) { int j = (i + 1) % loop.numVerts; - kadd(&k, cagnoli(loop.verts[i], loop.verts[j])); + kadd(&adder, cagnoli(loop.verts[i], loop.verts[j])); } - if (k.sum < 0.0) { - kadd(&k, 4.0 * M_PI); + if (adder.sum < 0.0) { + kadd(&adder, 4.0 * M_PI); } - *out = kresult(k); + *out = kresult(adder); return E_SUCCESS; } From 7528234979b9e5492da5bc1effe0199b9b4df595 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 22:25:10 -0500 Subject: [PATCH 15/74] bench --- src/apps/benchmarks/benchmarkArea.c | 7 ++----- src/h3lib/lib/area.c | 4 ++++ 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 69ce545cfb..87a05e1e0a 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,10 +9,6 @@ #include "iterators.h" #include "utility.h" -// TOOD: demonstrate the area alg works for global polygons. - -// TODO: is kadd faster without pointers? k = kadd(k, cellArea); - static void do_res_sum(int res) { Adder adder = {0, 0}; double cellArea; @@ -31,7 +27,8 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - for (int i = 0; i <= 6; i++) { + int MAX_RES = 8; + for (int i = 0; i <= MAX_RES; i++) { do_res_sum(i); } }); diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index a10e337dcf..6af6632122 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,6 +11,10 @@ #include "constants.h" #include "h3api.h" +// TOOD: demonstrate the area alg works for global polygons. + +// TODO: is kadd faster without pointers? k = kadd(k, cellArea); + static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; From 4a4fb073b4eebcd9237a0489e8fea1b717ccee3b Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 22:48:05 -0500 Subject: [PATCH 16/74] adder docs --- src/h3lib/include/adder.h | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index e391987964..9b8ad25e07 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -2,15 +2,16 @@ #define ADDER_H /* -This is a header-only implementation of a "compensated summation" algorithm, -which allows us to add up sequences of floating point numbers with better -precision than naive summation, especially when the terms in the sum -can vary in magnitude. +Header-only implementation of a "compensated summation" algorithm (Kahan +summation), which allows us to add up sequences of floating-point numbers +with better precision than naive summation, especially when the terms in +the sum vary significantly in magnitude. See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm -This is useful when computing the area of (multi-)polygons, which can involve -adding many small terms to large aggregates. -For example, D3 uses this when computing polygonal areas via the `Kahan` class: +This is useful when computing the area of (multi)polygons, which +often involves adding many small terms to a large aggregate. For example, +D3 uses an improved accuracy summation when computing polygonal areas via its +`Adder` class: https://github.com/d3/d3-geo/blob/main/src/area.js There are a few potential algorithms we might consider for summation: @@ -18,14 +19,14 @@ There are a few potential algorithms we might consider for summation: 1. Naive sum 2. Kahan summation 3. Neumaier summation -4. Other approaches like pairwise summation, or Pythons `fsum` +4. Other approaches like pairwise summation, or Python's `fsum` -We considered the first three for simplicity, and settled on Kahan summation, -as it acheives noticably better accuracy than naive summation and almost as good -accuracy as Neumaier, while only being slightly slower than naive, but slightly -faster and simpler than Neumaier. +We considered the first three for simplicity, and settled on Kahan summation: +it achieves noticeably better accuracy than naive summation and almost as +good accuracy as Neumaier, while being only slightly slower than naive and +slightly faster and simpler than Neumaier. -See also: https://github.com/python/cpython/issues/100425 for discussions of +See also: https://github.com/python/cpython/issues/100425 for discussion of tradeoffs between Kahan, Neumaier, and `fsum`. */ @@ -62,15 +63,15 @@ static inline double neumaier_result(Adder adder) { } static inline void kadd(Adder *adder, double x) { - // simple_add(adder, x); + simple_add(adder, x); // kahan_add(adder, x); - neumaier_add(adder, x); + // neumaier_add(adder, x); } static inline double kresult(Adder adder) { - // return simple_result(adder); + return simple_result(adder); // return kahan_result(adder); - return neumaier_result(adder); + // return neumaier_result(adder); } #endif From 6bd329904066292c54ee395de901b651c1d3fde3 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:03:18 -0500 Subject: [PATCH 17/74] Settled on Kahan implementation --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/apps/testapps/testH3CellAreaExhaustive.c | 2 +- src/h3lib/include/adder.h | 34 +------------------- src/h3lib/lib/area.c | 2 +- 4 files changed, 4 insertions(+), 36 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 87a05e1e0a..e9fa852d2d 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -19,7 +19,7 @@ static void do_res_sum(int res) { kadd(&adder, cellArea); } - double total_area = kresult(adder); + double total_area = adder.sum; double diff = fabs(total_area - 4 * M_PI); printf("res: %d, diff: %e\n", res, diff); } diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index b464014b46..4b0dc2b6b3 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -146,7 +146,7 @@ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), kadd(&adder, cellArea); } - t_assert(fabs(kresult(adder) - target) < tol, + t_assert(fabs(adder.sum - target) < tol, "sum of all cells should give earth area"); } diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 9b8ad25e07..1a6fc79695 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -35,43 +35,11 @@ typedef struct { double c; // compensation term } Adder; -static inline void simple_add(Adder *adder, double x) { adder->sum += x; } - -static inline void kahan_add(Adder *adder, double x) { +static inline void kadd(Adder *adder, double x) { double y = x - adder->c; double t = adder->sum + y; adder->c = (t - adder->sum) - y; adder->sum = t; } -static inline void neumaier_add(Adder *adder, double x) { - double t = adder->sum + x; - - if (fabs(adder->sum) >= fabs(x)) { - adder->c += (adder->sum - t) + x; - } else { - adder->c += (x - t) + adder->sum; - } - - adder->sum = t; -} - -static inline double simple_result(Adder adder) { return adder.sum; } -static inline double kahan_result(Adder adder) { return adder.sum; } -static inline double neumaier_result(Adder adder) { - return adder.sum + adder.c; -} - -static inline void kadd(Adder *adder, double x) { - simple_add(adder, x); - // kahan_add(adder, x); - // neumaier_add(adder, x); -} - -static inline double kresult(Adder adder) { - return simple_result(adder); - // return kahan_result(adder); - // return neumaier_result(adder); -} - #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 6af6632122..e026ddd525 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -41,7 +41,7 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { kadd(&adder, 4.0 * M_PI); } - *out = kresult(adder); + *out = adder.sum; return E_SUCCESS; } From 1b21bed993ef58a0ae24301fa91a8c0b6433c599 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:07:42 -0500 Subject: [PATCH 18/74] Tighten tolerances in testH3CellAreaExhaustive.c due to compensated sum --- src/apps/testapps/testH3CellAreaExhaustive.c | 26 +++++++++----------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 4b0dc2b6b3..315eb1e91e 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -178,28 +178,24 @@ SUITE(h3CellAreaExhaustive) { double km2 = rads2 * EARTH_RADIUS_KM * EARTH_RADIUS_KM; double m2 = km2 * 1000 * 1000; - // Notice the drop in accuracy at resolution 1. - // I think this has something to do with Class II vs Class III - // resolutions. - earth_area_test(0, H3_EXPORT(cellAreaRads2), rads2, 1e-14); earth_area_test(0, H3_EXPORT(cellAreaKm2), km2, 1e-6); earth_area_test(0, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(1, H3_EXPORT(cellAreaRads2), rads2, 1e-9); - earth_area_test(1, H3_EXPORT(cellAreaKm2), km2, 1e-1); - earth_area_test(1, H3_EXPORT(cellAreaM2), m2, 1e5); + earth_area_test(1, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(1, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(1, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(2, H3_EXPORT(cellAreaRads2), rads2, 1e-12); - earth_area_test(2, H3_EXPORT(cellAreaKm2), km2, 1e-5); + earth_area_test(2, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(2, H3_EXPORT(cellAreaKm2), km2, 1e-6); earth_area_test(2, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(3, H3_EXPORT(cellAreaRads2), rads2, 1e-11); - earth_area_test(3, H3_EXPORT(cellAreaKm2), km2, 1e-3); - earth_area_test(3, H3_EXPORT(cellAreaM2), m2, 1e3); + earth_area_test(3, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(3, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(3, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(4, H3_EXPORT(cellAreaRads2), rads2, 1e-11); - earth_area_test(4, H3_EXPORT(cellAreaKm2), km2, 1e-3); - earth_area_test(4, H3_EXPORT(cellAreaM2), m2, 1e2); + earth_area_test(4, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(4, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(4, H3_EXPORT(cellAreaM2), m2, 1e0); } } From 64853b3c011c6af823d5507cf6afc44406659513 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:31:33 -0500 Subject: [PATCH 19/74] 6 --- src/apps/benchmarks/benchmarkArea.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index e9fa852d2d..a29a771680 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -27,7 +27,7 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - int MAX_RES = 8; + int MAX_RES = 6; for (int i = 0; i <= MAX_RES; i++) { do_res_sum(i); } From 3547c47d829e2bf611d8fb4f7313f7f87b373a62 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:33:51 -0500 Subject: [PATCH 20/74] notes --- src/h3lib/include/adder.h | 6 ++++++ src/h3lib/lib/area.c | 2 -- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 1a6fc79695..53e31ac5f8 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -28,6 +28,12 @@ slightly faster and simpler than Neumaier. See also: https://github.com/python/cpython/issues/100425 for discussion of tradeoffs between Kahan, Neumaier, and `fsum`. + +## Usage + +Adder adder = {0.0, 0.0}; // initialize +kadd(&adder, x); // add x to summation +out = adder.sum; // extract final sum */ typedef struct { diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index e026ddd525..0a43bcf32a 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -13,8 +13,6 @@ // TOOD: demonstrate the area alg works for global polygons. -// TODO: is kadd faster without pointers? k = kadd(k, cellArea); - static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; From 0c73921821d631af1d159ab7fc21000db9348612 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 09:36:28 -0500 Subject: [PATCH 21/74] starting new tests --- CMakeLists.txt | 1 + CMakeTests.cmake | 1 + justfile | 3 ++- src/apps/testapps/testGeoLoopArea.c | 32 +++++++++++++++++++++++++++++ src/h3lib/lib/area.c | 2 ++ 5 files changed, 38 insertions(+), 1 deletion(-) create mode 100644 src/apps/testapps/testGeoLoopArea.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e4f88f7d9..dcc7069a46 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -280,6 +280,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testGeoLoopArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 1901a61fac..3508ff23a3 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,6 +259,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(testGeoLoopArea src/apps/testapps/testGeoLoopArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index 5d7ff258c8..8df37001ba 100644 --- a/justfile +++ b/justfile @@ -26,7 +26,8 @@ test: build # ./build/bin/testDirectedEdge # ./build/bin/testArea # just test-slow - just test-fast + # just test-fast + ./build/bin/testGeoLoopArea time: time ./build/bin/testArea diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c new file mode 100644 index 0000000000..f423c1db75 --- /dev/null +++ b/src/apps/testapps/testGeoLoopArea.c @@ -0,0 +1,32 @@ +/* + * Copyright 2020-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** @file + * @brief tests H3 cell area functions on a few specific cases + * + * usage: `testH3CellArea` + */ + +#include + +#include "h3Index.h" +#include "test.h" + +SUITE(h3CellArea) { + TEST(specific_cell_area) { + t_assert(false, "should fail"); + } +} diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 0a43bcf32a..bdbbc0f351 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -14,6 +14,7 @@ // TOOD: demonstrate the area alg works for global polygons. static inline double cagnoli(LatLng x, LatLng y) { + // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -74,5 +75,6 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { H3_MEMORY(free)(loop->verts); + loop->numVerts = 0; loop->verts = NULL; } From 31fd41644f3576efc1d22f97fdfc6ffc67a14627 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 18:49:08 -0500 Subject: [PATCH 22/74] tests work --- src/apps/testapps/testGeoLoopArea.c | 141 +++++++++++++++++++++++++++- 1 file changed, 138 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index f423c1db75..985a48b324 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -24,9 +24,144 @@ #include "h3Index.h" #include "test.h" +#include "utility.h" -SUITE(h3CellArea) { - TEST(specific_cell_area) { - t_assert(false, "should fail"); +// GeoMultiPolygon alloc_global_mpoly() { +// const int numPolygons = 8; +// const int numVerts = 3; +// const LatLng verts[numPolygons][numVerts] = { +// {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, +// {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, +// {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, +// {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, +// {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, +// {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, +// {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, +// {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, +// }; + +// GeoMultiPolygon mpoly = { +// .numPolygons = numPolygons, +// .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), +// }; + +// for (int i = 0; i < numPolygons; i++) { +// GeoPolygon *poly = &mpoly.polygons[i]; +// poly->numHoles = 0; +// poly->holes = NULL; +// poly->geoloop.numVerts = numVerts; +// poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); + +// for (int j = 0; j < numVerts; j++) { +// poly->geoloop.verts[j] = verts[i][j]; +// } +// } + +// qsort(mpoly.polygons, numPolygons, sizeof(GeoPolygon), cmp_poly_area); + +// return mpoly; +// } + +#define TOL 1e-14 + +static void helper(LatLng *verts, int numVerts, double target_area) { + GeoLoop loop = {.verts = verts, .numVerts = numVerts}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - target_area) < TOL, "area should match"); +} + +SUITE(geoLoopArea) { + TEST(triangle) { + // GeoLoop representing a triangle covering 1/8 of the globe, with + // points ordered according to right-hand rule (counter-clockwise). + // Expected area: pi/2 + LatLng verts[] = {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - M_PI / 2.0) < 1e-14, "area should match"); + } + TEST(reverse_triangle) { + // Reversed the order of the points, so they are in clockwise order. + // This GeoLoop represents the whole globe minus the triangle above. + // Expected area: 7*pi/2 + LatLng verts[] = {{0.0, M_PI_2}, {0.0, 0.0}, {M_PI_2, 0.0}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - 7 * M_PI / 2.0) < 1e-14, "area should match"); + } + + TEST(slice) { + // 1/4 slice of the globe, from north to south pole + // Expected area: pi + LatLng verts[] = {{M_PI_2, 0.0}, + {0.0, 0.0}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, + {M_PI_2, 0.0}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - M_PI) < 1e-14, "area should match"); + } + + TEST(reverse_slice) { + // 3/4 slice of the globe, from north to south pole, formed by reversing + // order of points from example above. + // Expected area: 3*pi + LatLng verts[] = {{M_PI_2, 0.0}, + {0.0, M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, 0.0}, + {M_PI_2, 0.0}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - 3 * M_PI) < 1e-14, "area should match"); + } + + TEST(hemisphere_east) { + // one hemisphere of globe. western or eastern? + // Expected area: 2*pi + LatLng verts[] = { + {M_PI_2, 0.0}, // north pole + {0.0, -M_PI_2}, // equator + {-M_PI_2, 0.0}, // south pole + {0.0, M_PI_2}, // equator + }; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - 2 * M_PI) < 1e-14, "area should match"); + } + + TEST(hemisphere_north) { + // one hemisphere of globe. northern + // Expected area: 2*pi + LatLng verts[] = { + {0.0, -M_PI}, // equator + {0.0, -M_PI_2}, // equator + {0.0, 0.0}, // equator + {0.0, M_PI_2}, // equator + // {0.0, M_PI}, // equator don't need last one here. last + // point + // unnecessary, but doesn't harm anything. + }; + + helper(verts, ARRAY_SIZE(verts), 2 * M_PI); } } From 61702f13dab6573b3a0dc2a18bc29041c2993db1 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 21:31:32 -0500 Subject: [PATCH 23/74] more tests --- src/apps/testapps/testGeoLoopArea.c | 194 +++++++++++++--------------- 1 file changed, 93 insertions(+), 101 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 985a48b324..a6663bd0a5 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -1,5 +1,5 @@ /* - * Copyright 2020-2021 Uber Technologies, Inc. + * Copyright 2025 Uber Technologies, Inc. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. @@ -15,9 +15,9 @@ */ /** @file - * @brief tests H3 cell area functions on a few specific cases + * @brief tests H3 GeoLoop area calculation * - * usage: `testH3CellArea` + * usage: `testGeoLoopArea` */ #include @@ -26,142 +26,134 @@ #include "test.h" #include "utility.h" -// GeoMultiPolygon alloc_global_mpoly() { -// const int numPolygons = 8; -// const int numVerts = 3; -// const LatLng verts[numPolygons][numVerts] = { -// {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, -// {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, -// {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, -// {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, -// {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, -// {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, -// {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, -// {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, -// }; - -// GeoMultiPolygon mpoly = { -// .numPolygons = numPolygons, -// .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), -// }; - -// for (int i = 0; i < numPolygons; i++) { -// GeoPolygon *poly = &mpoly.polygons[i]; -// poly->numHoles = 0; -// poly->holes = NULL; -// poly->geoloop.numVerts = numVerts; -// poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); - -// for (int j = 0; j < numVerts; j++) { -// poly->geoloop.verts[j] = verts[i][j]; -// } -// } - -// qsort(mpoly.polygons, numPolygons, sizeof(GeoPolygon), cmp_poly_area); - -// return mpoly; -// } - -#define TOL 1e-14 - -static void helper(LatLng *verts, int numVerts, double target_area) { +static void _compareArea(LatLng *verts, int numVerts, double target_area) { + double tol = 1e-14; GeoLoop loop = {.verts = verts, .numVerts = numVerts}; double out; t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - - t_assert(fabs(out - target_area) < TOL, "area should match"); + t_assert(fabs(out - target_area) < tol, "area should match"); } SUITE(geoLoopArea) { TEST(triangle) { // GeoLoop representing a triangle covering 1/8 of the globe, with // points ordered according to right-hand rule (counter-clockwise). - // Expected area: pi/2 - LatLng verts[] = {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {0.0, M_PI_2}, + }; - t_assert(fabs(out - M_PI / 2.0) < 1e-14, "area should match"); + _compareArea(verts, ARRAY_SIZE(verts), M_PI / 2); } + TEST(reverse_triangle) { // Reversed the order of the points, so they are in clockwise order. // This GeoLoop represents the whole globe minus the triangle above. - // Expected area: 7*pi/2 - LatLng verts[] = {{0.0, M_PI_2}, {0.0, 0.0}, {M_PI_2, 0.0}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + LatLng verts[] = { + {0.0, M_PI_2}, + {0.0, 0.0}, + {M_PI_2, 0.0}, + }; - t_assert(fabs(out - 7 * M_PI / 2.0) < 1e-14, "area should match"); + _compareArea(verts, ARRAY_SIZE(verts), 7 * M_PI / 2); } TEST(slice) { // 1/4 slice of the globe, from north to south pole - // Expected area: pi - LatLng verts[] = {{M_PI_2, 0.0}, - {0.0, 0.0}, - {-M_PI_2, 0.0}, - {0.0, M_PI_2}, - {M_PI_2, 0.0}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - - t_assert(fabs(out - M_PI) < 1e-14, "area should match"); + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), M_PI); } TEST(reverse_slice) { // 3/4 slice of the globe, from north to south pole, formed by reversing // order of points from example above. - // Expected area: 3*pi - LatLng verts[] = {{M_PI_2, 0.0}, - {0.0, M_PI_2}, - {-M_PI_2, 0.0}, - {0.0, 0.0}, - {M_PI_2, 0.0}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - - t_assert(fabs(out - 3 * M_PI) < 1e-14, "area should match"); + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, 0.0}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), 3 * M_PI); } TEST(hemisphere_east) { // one hemisphere of globe. western or eastern? - // Expected area: 2*pi LatLng verts[] = { - {M_PI_2, 0.0}, // north pole - {0.0, -M_PI_2}, // equator - {-M_PI_2, 0.0}, // south pole - {0.0, M_PI_2}, // equator + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, }; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - t_assert(fabs(out - 2 * M_PI) < 1e-14, "area should match"); + _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); } TEST(hemisphere_north) { // one hemisphere of globe. northern - // Expected area: 2*pi LatLng verts[] = { - {0.0, -M_PI}, // equator - {0.0, -M_PI_2}, // equator - {0.0, 0.0}, // equator - {0.0, M_PI_2}, // equator - // {0.0, M_PI}, // equator don't need last one here. last - // point - // unnecessary, but doesn't harm anything. + {0.0, -M_PI}, + {0.0, -M_PI_2}, + {0.0, 0.0}, + {0.0, M_PI_2}, }; - helper(verts, ARRAY_SIZE(verts), 2 * M_PI); + _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); + } + + TEST(percentage_slice) { + // Demonstrate that edge arcs between points in a polygon or geoloop + // should be less than 180 degrees (M_PI radians). + // + // Create a slice from north pole to equator and back to the north pole + // that sweeps out an edge arc of t * M_PI radians along the equator, + // so it should have an area of t*M_PI for t \in [0,1]. + // + // However, notice a discontinuity at t = 1 (i.e., M_PI radians or 180 + // degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. + // + // Recall that the area in stradians of the entire globe is 4*M_PI. + for (double t = 0; t <= 1.2; t += 0.01) { + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {0.0, t * M_PI - M_PI_2}, + }; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + double tol = 1e-13; + if (t < 0.99) { + // When t < 1, the largest angle in the triangle is less than + // 180 degrees + t_assert(fabs(out - t * M_PI) <= tol, "expected area"); + } else if (t > 1.01) { + // Discontinuity at t == 1. For t > 1, the triangle "flips", + // because the shortest geodesic path is on the other side of + // the globe. The triangle is now oriented in clockwise order, + // and the area computed is the area *outside* of the triangle, + // which starts at 3*pi. + t_assert(fabs(out - (2 + t) * M_PI) <= tol, "expected area"); + } + // Note that we avoid testing t == 1, since the triangle + // isn't well defined because there are many possible geodesic + // shortest paths when consecutive points are antipodal + // (180 degrees apart). + + // Also note that a large polygon with t > 1 is *still* + // represntable and we can compute its area accurately. + // We just need to add intermediate vertices so that + // no edge arc is greater than 180 degrees. + } } } From 4b33b0278a9cc0496275842f4e55176a7e0be5ac Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 21:57:57 -0500 Subject: [PATCH 24/74] comments --- src/apps/testapps/testGeoLoopArea.c | 110 ++++++++++++++++++---------- 1 file changed, 70 insertions(+), 40 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index a6663bd0a5..61b52d09c5 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -36,9 +36,18 @@ static void _compareArea(LatLng *verts, int numVerts, double target_area) { } SUITE(geoLoopArea) { - TEST(triangle) { - // GeoLoop representing a triangle covering 1/8 of the globe, with - // points ordered according to right-hand rule (counter-clockwise). + TEST(triangle_basic) { + /* + GeoLoop representing a triangle covering 1/8 of the globe, with + points ordered according to right-hand rule (counter-clockwise). + + The triangle starts at the north pole, moves down 90 degrees to the + equator, and then sweeps out 90 degrees along the equator + before returning to the north pole. + + The globe has an area of 4*pi radians, so this 1/8 triangle piece of + the globe should have area pi/2. + */ LatLng verts[] = { {M_PI_2, 0.0}, {0.0, 0.0}, @@ -48,9 +57,14 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), M_PI / 2); } - TEST(reverse_triangle) { - // Reversed the order of the points, so they are in clockwise order. - // This GeoLoop represents the whole globe minus the triangle above. + TEST(triangle_reversed) { + /* + Reverse the order of the points in the triangle from the previous test, + so that they are in clockwise order. + + Since the points are in clockwise order, GeoLoop represents the whole + globe minus the triangle above. + */ LatLng verts[] = { {0.0, M_PI_2}, {0.0, 0.0}, @@ -61,7 +75,11 @@ SUITE(geoLoopArea) { } TEST(slice) { - // 1/4 slice of the globe, from north to south pole + /* + Stitch two 1/8 triangles together, sharing an edge along the equator + to create a 1/4 slice of the globe, with vertices at the north and south + pole. + */ LatLng verts[] = { {M_PI_2, 0.0}, {0.0, 0.0}, @@ -72,9 +90,11 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), M_PI); } - TEST(reverse_slice) { - // 3/4 slice of the globe, from north to south pole, formed by reversing - // order of points from example above. + TEST(slice_reversed) { + /* + 3/4 slice of the globe, from north to south pole, formed by reversing + order of points from example above. + */ LatLng verts[] = { {M_PI_2, 0.0}, {0.0, M_PI_2}, @@ -86,19 +106,23 @@ SUITE(geoLoopArea) { } TEST(hemisphere_east) { - // one hemisphere of globe. western or eastern? + /* + Stitch 4 1/8 triangles together to cover the eastern hemisphere. + */ LatLng verts[] = { {M_PI_2, 0.0}, - {0.0, -M_PI_2}, + {0.0, 0}, {-M_PI_2, 0.0}, - {0.0, M_PI_2}, + {0.0, M_PI}, }; _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); } TEST(hemisphere_north) { - // one hemisphere of globe. northern + /* + Stitch 4 1/8 triangles together to cover the northern hemisphere. + */ LatLng verts[] = { {0.0, -M_PI}, {0.0, -M_PI_2}, @@ -109,18 +133,20 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); } - TEST(percentage_slice) { - // Demonstrate that edge arcs between points in a polygon or geoloop - // should be less than 180 degrees (M_PI radians). - // - // Create a slice from north pole to equator and back to the north pole - // that sweeps out an edge arc of t * M_PI radians along the equator, - // so it should have an area of t*M_PI for t \in [0,1]. - // - // However, notice a discontinuity at t = 1 (i.e., M_PI radians or 180 - // degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. - // - // Recall that the area in stradians of the entire globe is 4*M_PI. + TEST(percentageSlice) { + /* + Demonstrate that edge arcs between points in a polygon or geoloop + should be less than 180 degrees (M_PI radians). + + Create a triangle from north pole to equator and back to the north pole + that sweeps out an edge arc of t * M_PI radians along the equator, + so it should have an area of t*M_PI for t \in [0,1]. + + However, there is a discontinuity at t = 1 (i.e., M_PI radians or 180 + degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. + + Recall that the area in stradians of the entire globe is 4*M_PI. + */ for (double t = 0; t <= 1.2; t += 0.01) { LatLng verts[] = { {M_PI_2, 0.0}, @@ -138,22 +164,26 @@ SUITE(geoLoopArea) { // 180 degrees t_assert(fabs(out - t * M_PI) <= tol, "expected area"); } else if (t > 1.01) { - // Discontinuity at t == 1. For t > 1, the triangle "flips", - // because the shortest geodesic path is on the other side of - // the globe. The triangle is now oriented in clockwise order, - // and the area computed is the area *outside* of the triangle, - // which starts at 3*pi. + /* + Discontinuity at t == 1. For t > 1, the triangle "flips", + because the shortest geodesic path is on the other side of + the globe. The triangle is now oriented in clockwise order, + and the area computed is the area *outside* of the triangle, + which starts at 3*pi. + */ t_assert(fabs(out - (2 + t) * M_PI) <= tol, "expected area"); } - // Note that we avoid testing t == 1, since the triangle - // isn't well defined because there are many possible geodesic - // shortest paths when consecutive points are antipodal - // (180 degrees apart). - - // Also note that a large polygon with t > 1 is *still* - // represntable and we can compute its area accurately. - // We just need to add intermediate vertices so that - // no edge arc is greater than 180 degrees. + /* + Note that we avoid testing t == 1, since the triangle + isn't well defined because there are many possible geodesic + shortest paths when consecutive points are antipodal + (180 degrees apart). + + Also note that a large polygon with t > 1 is *still* + representable and we can compute its area accurately; + we just need to add intermediate vertices so that + no edge arc is greater than 180 degrees. + */ } } } From 26eb14376f1916f65d9996774168b66d1f732798 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 22:45:06 -0500 Subject: [PATCH 25/74] docstring --- src/h3lib/lib/area.c | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index bdbbc0f351..a1c00a8526 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,8 +11,6 @@ #include "constants.h" #include "h3api.h" -// TOOD: demonstrate the area alg works for global polygons. - static inline double cagnoli(LatLng x, LatLng y) { // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 x.lat = x.lat / 2.0 + M_PI / 4.0; @@ -28,6 +26,30 @@ static inline double cagnoli(LatLng x, LatLng y) { return -2.0 * atan2(sa * sd, sa * cd + ca); } +/** + * Area in radians^2 enclosed by vertices in GeoLoop. + * + * The GeoLoop should represent a simple curve with no self-intersections. + * Vertices should be ordered according to the "right hand rule". + * That is, if you are looking from outer space at a spherical polygon + * on the surface of the earth whose interior is contained within a hemisphere, + * then the vertices should be ordered counter-clockwise. The interior of the + * polygon is to the left of a person walking along the boundary of the polygon + * in the counter-clockwise direction. + * + * Note that GeoLoops do not need to repeat the first vertex at the end of the + * array to close the loop; this is done automatically. + * + * The area of the entire globe is 4*pi radians^2. If, for example, you have a + * small GeoLoop with area `a << 4*pi` and then reverse the order of the + * vertices, you produce GeoLoop with area `4*pi - a`, since by the right hand + * rule, the interior of the new loop's interior is the majority of the globe, + * or "everything except the original polygon". + * + * @param loop GeoLoop of boundary vertices in counter-clockwise order + * @param out loop area in radians^2, in interval [0, 4*pi] + * @return E_SUCCESS on success, or an error code otherwise + */ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { Adder adder = {0.0, 0.0}; @@ -47,11 +69,7 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { /** * Area of H3 cell in radians^2. * - * The area is calculated by breaking the cell into spherical triangles and - * summing up their areas. Note that some H3 cells (hexagons and pentagons) - * are irregular, and have more than 6 or 5 sides. - * - * todo: optimize the computation by re-using the edges shared between triangles + * Uses `geoLoopArea` to compute cell area. * * @param cell H3 cell * @param out cell area in radians^2 @@ -72,9 +90,3 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { return E_SUCCESS; } - -void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { - H3_MEMORY(free)(loop->verts); - loop->numVerts = 0; - loop->verts = NULL; -} From 06b2a6b8b6d9348311483d296b5af740d584623a Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:08:53 -0500 Subject: [PATCH 26/74] better docstrings --- src/h3lib/lib/area.c | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index a1c00a8526..db640bc8a8 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,8 +11,11 @@ #include "constants.h" #include "h3api.h" + +// Cagnoli contribution for edge arc x to y, following d3-geo’s +// area implementation: +// https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 static inline double cagnoli(LatLng x, LatLng y) { - // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -31,20 +34,30 @@ static inline double cagnoli(LatLng x, LatLng y) { * * The GeoLoop should represent a simple curve with no self-intersections. * Vertices should be ordered according to the "right hand rule". - * That is, if you are looking from outer space at a spherical polygon + * That is, if you are looking from outer space at a spherical polygon loop * on the surface of the earth whose interior is contained within a hemisphere, * then the vertices should be ordered counter-clockwise. The interior of the - * polygon is to the left of a person walking along the boundary of the polygon + * loop is to the left of a person walking along the boundary of the polygon * in the counter-clockwise direction. * * Note that GeoLoops do not need to repeat the first vertex at the end of the * array to close the loop; this is done automatically. * + * The edge arcs between adjacent vertices are assumed to be the shortest + * geodesic path between them; that is, all arcs are interpreted to be less + * than 180 degrees or pi radians. + * Avoid arcs that are exactly pi (i.e, two antipodal vertices). + * "Large" polygon loops (e.g., that cannot be contained in a hemisphere) can + * still be constructed by using intermediate vertices with arcs less than + * 180 degrees, and the loop area will still be computed correctly. + * * The area of the entire globe is 4*pi radians^2. If, for example, you have a * small GeoLoop with area `a << 4*pi` and then reverse the order of the - * vertices, you produce GeoLoop with area `4*pi - a`, since by the right hand - * rule, the interior of the new loop's interior is the majority of the globe, + * vertices, you produce GeoLoop with area `4*pi - a`, since, by the right hand + * rule, the new loop's interior is the majority of the globe, * or "everything except the original polygon". + * Note that the area enclosed by the loop is determined by the vertex order; + * this function does **not** return `min(a, 4*pi - a)`. * * @param loop GeoLoop of boundary vertices in counter-clockwise order * @param out loop area in radians^2, in interval [0, 4*pi] @@ -58,6 +71,10 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { kadd(&adder, cagnoli(loop.verts[i], loop.verts[j])); } + // The Cagnoli sum above yields a signed area, with the sign switching + // with the orientation of the vertices. Since we want our area to always be + // positive, we normalize into [0, 4*pi] by adding 4*pi when the signed + // area is negative. if (adder.sum < 0.0) { kadd(&adder, 4.0 * M_PI); } From f55e1b14366a6fbc085e3cffb6f20f2b68fb34e1 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:10:54 -0500 Subject: [PATCH 27/74] Adder note --- src/h3lib/lib/area.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index db640bc8a8..ac8f82a3d4 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,7 +11,6 @@ #include "constants.h" #include "h3api.h" - // Cagnoli contribution for edge arc x to y, following d3-geo’s // area implementation: // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 @@ -64,6 +63,8 @@ static inline double cagnoli(LatLng x, LatLng y) { * @return E_SUCCESS on success, or an error code otherwise */ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { + // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli + // terms Adder adder = {0.0, 0.0}; for (int i = 0; i < loop.numVerts; i++) { From 0c939fce58de248d6424c1197effcb3a8d636c03 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:20:11 -0500 Subject: [PATCH 28/74] don't need these imports --- justfile | 4 ++-- src/apps/testapps/testGeoLoopArea.c | 2 +- src/apps/testapps/testH3CellArea.c | 2 +- src/h3lib/lib/area.c | 7 ------- 4 files changed, 4 insertions(+), 11 deletions(-) diff --git a/justfile b/justfile index 8df37001ba..a9919cc8ce 100644 --- a/justfile +++ b/justfile @@ -26,8 +26,8 @@ test: build # ./build/bin/testDirectedEdge # ./build/bin/testArea # just test-slow - # just test-fast - ./build/bin/testGeoLoopArea + just test-fast + # ./build/bin/testGeoLoopArea time: time ./build/bin/testArea diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 61b52d09c5..ffac0a505f 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -22,7 +22,7 @@ #include -#include "h3Index.h" +#include "h3api.h" #include "test.h" #include "utility.h" diff --git a/src/apps/testapps/testH3CellArea.c b/src/apps/testapps/testH3CellArea.c index 060e14d169..4189dfc7de 100644 --- a/src/apps/testapps/testH3CellArea.c +++ b/src/apps/testapps/testH3CellArea.c @@ -22,7 +22,7 @@ #include -#include "h3Index.h" +#include "h3api.h" #include "test.h" static const double areasKm2[] = { diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index ac8f82a3d4..ef1bb845c1 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -1,13 +1,6 @@ -#include #include -#include -#include -#include -#include -#include #include "adder.h" -#include "alloc.h" #include "constants.h" #include "h3api.h" From a292924a53d230ac6351a08c5997421292dfcbfd Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:24:34 -0500 Subject: [PATCH 29/74] clean a few more imports --- src/apps/benchmarks/benchmarkArea.c | 3 --- src/h3lib/lib/area.c | 19 +++++++++++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index a29a771680..c353ed4628 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -1,5 +1,3 @@ -#include -#include #include #include @@ -7,7 +5,6 @@ #include "benchmark.h" #include "h3api.h" #include "iterators.h" -#include "utility.h" static void do_res_sum(int res) { Adder adder = {0, 0}; diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index ef1bb845c1..73c5e819ac 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -1,3 +1,22 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file area.c + * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, cells, polygons, multipolygons, etc.) + */ + #include #include "adder.h" From 00f45a2f0f8ba7628d0d09853389caf204b6f188 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:24:55 -0500 Subject: [PATCH 30/74] format --- src/h3lib/lib/area.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 73c5e819ac..88cec9565d 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -14,7 +14,8 @@ * limitations under the License. */ /** @file area.c - * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, cells, polygons, multipolygons, etc.) + * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, + * cells, polygons, multipolygons, etc.) */ #include From 39be84c2c6c32e4f307ca62ed95cf8f12aadc732 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:30:05 -0500 Subject: [PATCH 31/74] needs constants --- src/apps/benchmarks/benchmarkArea.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index c353ed4628..997e95e69c 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -5,6 +5,7 @@ #include "benchmark.h" #include "h3api.h" #include "iterators.h" +#include "constants.h" static void do_res_sum(int res) { Adder adder = {0, 0}; From 53e5be9ccbd3b7102699618f088be20cd3d3b646 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:35:04 -0500 Subject: [PATCH 32/74] format --- src/apps/benchmarks/benchmarkArea.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 997e95e69c..c11f4842b1 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -3,9 +3,9 @@ #include "adder.h" #include "benchmark.h" +#include "constants.h" #include "h3api.h" #include "iterators.h" -#include "constants.h" static void do_res_sum(int res) { Adder adder = {0, 0}; From 139f483147413e0b1affcc60e642ec94d73b4893 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 09:53:58 -0500 Subject: [PATCH 33/74] try fuzzer --- CMakeLists.txt | 2 ++ src/apps/fuzzers/README.md | 1 + src/apps/fuzzers/fuzzerGeoLoopArea.c | 39 ++++++++++++++++++++++++++++ src/h3lib/include/h3api.h.in | 1 - 4 files changed, 42 insertions(+), 1 deletion(-) create mode 100644 src/apps/fuzzers/fuzzerGeoLoopArea.c diff --git a/CMakeLists.txt b/CMakeLists.txt index dcc7069a46..d2bea788f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -310,6 +310,7 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerCellToChildPos.c src/apps/fuzzers/fuzzerInternalAlgos.c src/apps/fuzzers/fuzzerInternalCoordIjk.c + src/apps/fuzzers/fuzzerGeoLoopArea.c src/apps/benchmarks/benchmarkPolygonToCells.c src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c src/apps/benchmarks/benchmarkPolygon.c @@ -639,6 +640,7 @@ if(BUILD_FUZZERS) add_h3_fuzzer(fuzzerPolygonToCellsExperimentalNoHoles src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c) add_h3_fuzzer(fuzzerCellToChildPos src/apps/fuzzers/fuzzerCellToChildPos.c) + add_h3_fuzzer(fuzzerGeoLoopArea src/apps/fuzzers/fuzzerGeoLoopArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_fuzzer(fuzzerInternalAlgos src/apps/fuzzers/fuzzerInternalAlgos.c) diff --git a/src/apps/fuzzers/README.md b/src/apps/fuzzers/README.md index 31135f0519..5b2a931d9f 100644 --- a/src/apps/fuzzers/README.md +++ b/src/apps/fuzzers/README.md @@ -64,6 +64,7 @@ The public API of H3 is covered in the following fuzzers: | originToDirectedEdges | [fuzzerDirectedEdge](./fuzzerDirectedEdge.c) | polygonToCells | [fuzzerPoylgonToCells](./fuzzerPolygonToCells.c) | polygonToCellsExperimental | [fuzzerPoylgonToCellsExperimental](./fuzzerPolygonToCellsExperimental.c) [fuzzerPoylgonToCellsExperimentalNoHoles](./fuzzerPolygonToCellsExperimentalNoHoles.c) +| geoLoopArea | [fuzzerGeoLoopArea](./fuzzerGeoLoopArea.c) | radsToDegs | Trivial | stringToH3 | [fuzzerIndexIO](./fuzzerIndexIO.c) | uncompactCells | [fuzzerCompact](./fuzzerCompact.c) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c new file mode 100644 index 0000000000..a5510d1519 --- /dev/null +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -0,0 +1,39 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file + * @brief Fuzzer program for geoLoopArea + */ + +#include "aflHarness.h" +#include "h3api.h" +#include "utility.h" + +int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { + // Interpret the input buffer as an array of LatLng vertices. + int numVerts = (int)(size / sizeof(LatLng)); + // Need at least 3 vertices to form a polygon. + if (numVerts < 3) { + return 0; + } + + GeoLoop loop = {.numVerts = numVerts, .verts = (LatLng *)data}; + double area; + H3_EXPORT(geoLoopArea)(loop, &area); + + return 0; +} + +AFL_HARNESS_MAIN(sizeof(LatLng) * 1024); diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 042b778471..8fe85cbcb0 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,7 +348,6 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop); DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); /** @defgroup degsToRads degsToRads From afa64f166931325f9112c614fb24f62e05144df0 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 10:24:08 -0500 Subject: [PATCH 34/74] .mdx --- website/docs/api/regions.mdx | 90 ++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 048d62935e..0d8ede3871 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -647,4 +647,94 @@ This function exists for memory management and is not exposed. +## geoLoopArea + +Compute the spherical area, in square radians, that is enclosed by a closed +`GeoLoop`. The loop must be a simple polygon traced counter-clockwise as viewed +from above the Earth, and each edge is interpreted as the shortest geodesic +between consecutive vertices (i.e., no side spans more than π radians). Areas +are normalized to the [0, 4π] interval, so reversing the vertex order yields +the complementary area of the loop. + + + + +```c +H3Error geoLoopArea(GeoLoop loop, double *out); +``` + +The caller supplies the `GeoLoop` memory (for example, by filling `loop.verts` +with previously allocated `LatLng` vertices) and receives the enclosed area in +square radians in `out`. Returns 0 (`E_SUCCESS`) on success. + + + + +:::note + +This function is not exposed in the Java bindings. + +::: + + + + +:::note + +This function is not exposed in the JavaScript bindings. + +::: + + + + +:::note + +This function is not exposed in the Python bindings. + +TODO: direct users to create a polygon with no holes. + +::: + + + + +:::note + +This function is not exposed in the Go bindings. + +::: + + + + +:::note + +This function is not exposed in the DuckDB bindings. + +::: + + + + +:::note + +This function is not exposed in the CLI bindings. + +::: + + + From f855776e60da80714c673e9a316d29e4697d966f Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 10:27:30 -0500 Subject: [PATCH 35/74] docs in h3api.h.in --- src/h3lib/include/h3api.h.in | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 8fe85cbcb0..19837bcb42 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,7 +348,21 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ +/** @defgroup geoLoopArea geoLoopArea + * Functions for geoLoopArea + * @{ + */ +/** + * @brief Area in radians^2 enclosed by vertices in a GeoLoop. + * + * Vertices should define a simple polygon ordered counter-clockwise according + * to the right-hand rule. Edge arcs follow the shortest geodesic path (less + * than pi radians), and the first vertex does not need to be repeated to close + * the loop. The returned area depends on vertex orientation and is in the + * interval [0, 4*pi]. + */ DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); +/** @} */ /** @defgroup degsToRads degsToRads * Functions for degsToRads From e48b735105053541eff39352cdc3c08a9e9a6d35 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:05:37 -0500 Subject: [PATCH 36/74] rads2 --- src/apps/fuzzers/fuzzerGeoLoopArea.c | 4 ++-- src/apps/testapps/testGeoLoopArea.c | 4 ++-- src/h3lib/include/h3api.h.in | 6 +++--- src/h3lib/lib/area.c | 6 +++--- website/docs/api/regions.mdx | 4 ++-- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index a5510d1519..88b96e1b7b 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -14,7 +14,7 @@ * limitations under the License. */ /** @file - * @brief Fuzzer program for geoLoopArea + * @brief Fuzzer program for geoLoopAreaRads2 */ #include "aflHarness.h" @@ -31,7 +31,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { GeoLoop loop = {.numVerts = numVerts, .verts = (LatLng *)data}; double area; - H3_EXPORT(geoLoopArea)(loop, &area); + H3_EXPORT(geoLoopAreaRads2)(loop, &area); return 0; } diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index ffac0a505f..f043cc3e7d 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -31,7 +31,7 @@ static void _compareArea(LatLng *verts, int numVerts, double target_area) { GeoLoop loop = {.verts = verts, .numVerts = numVerts}; double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); t_assert(fabs(out - target_area) < tol, "area should match"); } @@ -156,7 +156,7 @@ SUITE(geoLoopArea) { GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); double tol = 1e-13; if (t < 0.99) { diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 19837bcb42..1220414f08 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,8 +348,8 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -/** @defgroup geoLoopArea geoLoopArea - * Functions for geoLoopArea +/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 + * Functions for geoLoopAreaRads2 * @{ */ /** @@ -361,7 +361,7 @@ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); * the loop. The returned area depends on vertex orientation and is in the * interval [0, 4*pi]. */ -DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); +DECLSPEC H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out); /** @} */ /** @defgroup degsToRads degsToRads diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 88cec9565d..f94b68f306 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -75,7 +75,7 @@ static inline double cagnoli(LatLng x, LatLng y) { * @param out loop area in radians^2, in interval [0, 4*pi] * @return E_SUCCESS on success, or an error code otherwise */ -H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { +H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms Adder adder = {0.0, 0.0}; @@ -100,7 +100,7 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { /** * Area of H3 cell in radians^2. * - * Uses `geoLoopArea` to compute cell area. + * Uses `geoLoopAreaRads2` to compute cell area. * * @param cell H3 cell * @param out cell area in radians^2 @@ -114,7 +114,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { } GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; - err = H3_EXPORT(geoLoopArea)(loop, out); + err = H3_EXPORT(geoLoopAreaRads2)(loop, out); if (err) { return err; } diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 0d8ede3871..7ccf0c1248 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -647,7 +647,7 @@ This function exists for memory management and is not exposed. -## geoLoopArea +## geoLoopAreaRads2 Compute the spherical area, in square radians, that is enclosed by a closed `GeoLoop`. The loop must be a simple polygon traced counter-clockwise as viewed @@ -672,7 +672,7 @@ the complementary area of the loop. ```c -H3Error geoLoopArea(GeoLoop loop, double *out); +H3Error geoLoopAreaRads2(GeoLoop loop, double *out); ``` The caller supplies the `GeoLoop` memory (for example, by filling `loop.verts` From 11b277a2b97282d52393d7c40bbdd40b468d94fb Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:25:47 -0500 Subject: [PATCH 37/74] python plan --- website/docs/api/regions.mdx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 7ccf0c1248..57711dfcde 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -702,9 +702,13 @@ This function is not exposed in the JavaScript bindings. :::note -This function is not exposed in the Python bindings. +This function isn't exposed directly, but can be used +via `h3.h3shape_area()`: -TODO: direct users to create a polygon with no holes. +```python +poly = h3.LatLngPoly([(90,0),(0,0),(0,90)]) +h3.h3shape_area(poly, unit='km^2') +``` ::: From 226da8cadffb27f0110e2d8dcf43d537d03e04f1 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:32:43 -0500 Subject: [PATCH 38/74] show, don't tell --- src/apps/testapps/testGeoLoopArea.c | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index f043cc3e7d..95c3c5f4aa 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -178,12 +178,25 @@ SUITE(geoLoopArea) { isn't well defined because there are many possible geodesic shortest paths when consecutive points are antipodal (180 degrees apart). - - Also note that a large polygon with t > 1 is *still* - representable and we can compute its area accurately; - we just need to add intermediate vertices so that - no edge arc is greater than 180 degrees. */ } } + + TEST(percentageSlice_large) { + /* + Continuing from the test above, note that a large polygon with + t > 1 is *still* representable and we can compute its area accurately; + we just need to add intermediate vertices so that + no edge arc is greater than 180 degrees. + */ + double t = 1.2; + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {0.0, 0.0}, // Extra vertex so every angle is < 180 degrees + {0.0, t * M_PI - M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), t * M_PI); + } } From c274bff131a6bb129194e2d75610e5d8d460e5e9 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:55:40 -0500 Subject: [PATCH 39/74] benchmark clean up --- src/apps/benchmarks/benchmarkArea.c | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index c11f4842b1..cf412ce3b5 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -1,4 +1,5 @@ #include +#include #include #include "adder.h" @@ -7,7 +8,7 @@ #include "h3api.h" #include "iterators.h" -static void do_res_sum(int res) { +static void doResSum(int res, bool print) { Adder adder = {0, 0}; double cellArea; IterCellsResolution iter = iterInitRes(res); @@ -19,15 +20,24 @@ static void do_res_sum(int res) { double total_area = adder.sum; double diff = fabs(total_area - 4 * M_PI); - printf("res: %d, diff: %e\n", res, diff); + if (print) { + printf("res: %d, diff: %e\n", res, diff); + } } BEGIN_BENCHMARKS(); -BENCHMARK(directedEdgeToBoundary, 1, { - int MAX_RES = 6; +static int MAX_RES = 4; + +BENCHMARK(allCellsAtRes_print, 1, { + for (int i = 0; i <= MAX_RES; i++) { + doResSum(i, true); + } +}); + +BENCHMARK(allCellsAtRes_noprint, 100, { for (int i = 0; i <= MAX_RES; i++) { - do_res_sum(i); + doResSum(i, false); } }); From 43f19a949bc1728dec3e545d5240d33a075d0746 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:04:19 -0500 Subject: [PATCH 40/74] never say never again --- src/h3lib/lib/area.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index f94b68f306..a62f5b68bf 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -115,7 +115,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; err = H3_EXPORT(geoLoopAreaRads2)(loop, out); - if (err) { + if (NEVER(err)) { return err; } From fffd72399e9956da2cd117fa0fe8eb997bff8c73 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:07:25 -0500 Subject: [PATCH 41/74] ugh, would have been cooler if it compiled the first time --- src/h3lib/lib/area.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index a62f5b68bf..9ef3b3c629 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -22,6 +22,7 @@ #include "adder.h" #include "constants.h" +#include "h3Assert.h" #include "h3api.h" // Cagnoli contribution for edge arc x to y, following d3-geo’s From efcbea0dae9b3b11f4abdf32a5c5e293ac479e5d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:25:51 -0500 Subject: [PATCH 42/74] remove justfile --- justfile | 48 ------------------------------------------------ 1 file changed, 48 deletions(-) delete mode 100644 justfile diff --git a/justfile b/justfile deleted file mode 100644 index a9919cc8ce..0000000000 --- a/justfile +++ /dev/null @@ -1,48 +0,0 @@ -# TODO: Remove before landing - -init: purge - mkdir build - -build: - cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - -profile: init - cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make - # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - xcrun xctrace record \ - --template 'Time Profiler' \ - --output h3-prof.trace \ - --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon - open h3-prof.trace - -purge: - rm -rf build - rm -rf *.trace - rm -rf .ipynb_checkpoints - -test: build - # ./build/bin/testH3CellAreaExhaustive - # ./build/bin/testEdgeCellsToPoly - # ./build/bin/testDirectedEdge - # ./build/bin/testArea - # just test-slow - just test-fast - # ./build/bin/testGeoLoopArea - -time: - time ./build/bin/testArea - -test-fast: build - cd build; make test-fast - -test-slow: build - cd build; make test - -bench: build - ./build/bin/benchmarkArea - # ./build/bin/benchmarkCellsToLinkedMultiPolygon - # ./build/bin/benchmarkCellsToMultiPolygon - # ./build/bin/benchmarkDirectedEdge - -fuzz: - cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers From 6e955c9541a0fbc45743fa6a1c50acb3a90c12a8 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:44:00 -0500 Subject: [PATCH 43/74] lighten the h3api.h.in description, as per usual in that file --- CHANGELOG.md | 5 +++++ src/h3lib/include/h3api.h.in | 8 +------- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 415bcdf886..a08417da08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ The public API of this library consists of the functions declared in file [h3api.h.in](./src/h3lib/include/h3api.h.in). ## [Unreleased] +### Added +- `geoLoopArea` function (#1101) + +### Changed +- `cellAreaRads2` uses `geoLoopArea` (#1101) ## [4.4.1] - 2025-11-11 ### Fixed diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 1220414f08..de589d50fa 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -353,13 +353,7 @@ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); * @{ */ /** - * @brief Area in radians^2 enclosed by vertices in a GeoLoop. - * - * Vertices should define a simple polygon ordered counter-clockwise according - * to the right-hand rule. Edge arcs follow the shortest geodesic path (less - * than pi radians), and the first vertex does not need to be repeated to close - * the loop. The returned area depends on vertex orientation and is in the - * interval [0, 4*pi]. + * @brief Area in radians^2 enclosed by vertices of a GeoLoop. */ DECLSPEC H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out); /** @} */ From dfeb6d828ca7963f7fd3cd08082c5e6c138f065e Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 18:23:08 -0500 Subject: [PATCH 44/74] simplify adder initialization --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/apps/testapps/testH3CellAreaExhaustive.c | 2 +- src/h3lib/include/adder.h | 2 +- src/h3lib/lib/area.c | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index cf412ce3b5..fe1ffa4e29 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,7 +9,7 @@ #include "iterators.h" static void doResSum(int res, bool print) { - Adder adder = {0, 0}; + Adder adder = {0}; double cellArea; IterCellsResolution iter = iterInitRes(res); diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 315eb1e91e..5a8bd30aa6 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -137,7 +137,7 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - Adder adder = {0.0, 0.0}; + Adder adder = {0}; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 53e31ac5f8..57682d4a78 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -31,7 +31,7 @@ tradeoffs between Kahan, Neumaier, and `fsum`. ## Usage -Adder adder = {0.0, 0.0}; // initialize +Adder adder = {0}; // initialize kadd(&adder, x); // add x to summation out = adder.sum; // extract final sum */ diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 9ef3b3c629..5b0b2ac3cb 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -79,7 +79,7 @@ static inline double cagnoli(LatLng x, LatLng y) { H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms - Adder adder = {0.0, 0.0}; + Adder adder = {0}; for (int i = 0; i < loop.numVerts; i++) { int j = (i + 1) % loop.numVerts; From 6db33a30962f327c1b94307f3869ab703623790f Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 20:33:07 -0500 Subject: [PATCH 45/74] drop numVerts < 3 --- justfile | 53 ++++++++++++++++++++++++++++ src/apps/fuzzers/fuzzerGeoLoopArea.c | 11 +++--- 2 files changed, 57 insertions(+), 7 deletions(-) create mode 100644 justfile diff --git a/justfile b/justfile new file mode 100644 index 0000000000..cffb5a6fde --- /dev/null +++ b/justfile @@ -0,0 +1,53 @@ +# TODO: Remove before landing + +init: purge + mkdir build + +build: + cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + +profile: init + cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make + # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + xcrun xctrace record \ + --template 'Time Profiler' \ + --output h3-prof.trace \ + --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon + open h3-prof.trace + +purge: + rm -rf build + rm -rf *.trace + rm -rf .ipynb_checkpoints + +test: build + # ./build/bin/testH3CellAreaExhaustive + # ./build/bin/testEdgeCellsToPoly + # ./build/bin/testDirectedEdge + # ./build/bin/testArea + # just test-slow + just test-fast + # ./build/bin/testGeoLoopArea + +time: + time ./build/bin/testArea + +test-fast: build + cd build; make test-fast + +test-slow: build + cd build; make test + +bench: build + ./build/bin/benchmarkArea + # ./build/bin/benchmarkCellsToLinkedMultiPolygon + # ./build/bin/benchmarkCellsToMultiPolygon + # ./build/bin/benchmarkDirectedEdge + +fuzz: + cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers + +coverage: purge + mkdir build + cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage + open build/coverage/index.html diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index 88b96e1b7b..b13b69f488 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -23,13 +23,10 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { // Interpret the input buffer as an array of LatLng vertices. - int numVerts = (int)(size / sizeof(LatLng)); - // Need at least 3 vertices to form a polygon. - if (numVerts < 3) { - return 0; - } - - GeoLoop loop = {.numVerts = numVerts, .verts = (LatLng *)data}; + GeoLoop loop = { + .numVerts = (int)(size / sizeof(LatLng)), + .verts = (LatLng *)data, + }; double area; H3_EXPORT(geoLoopAreaRads2)(loop, &area); From e2c20b0d96e2c96997be7be822cfb713cff6cf19 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 20:42:05 -0500 Subject: [PATCH 46/74] degenerate loop tests --- src/apps/testapps/testGeoLoopArea.c | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 95c3c5f4aa..5de08002a5 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -199,4 +199,27 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), t * M_PI); } + + TEST(degenerateLoop2) { + // TODO: should we raise an error when given degenerate loops? + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + }; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } + + TEST(degenerateLoop1) { + // TODO: should we raise an error when given degenerate loops? + LatLng verts[] = { + {0.0, 0.0}, + }; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } + + TEST(degenerateLoop0) { + // TODO: should we raise an error when given degenerate loops? + LatLng verts[] = {}; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } } From eeb3046fc694ec2bd310edf225fac47343b630df Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 20:51:44 -0500 Subject: [PATCH 47/74] try _compareArea(NULL, 0, 0.0); --- src/apps/testapps/testGeoLoopArea.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 5de08002a5..db331c71f8 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -219,7 +219,6 @@ SUITE(geoLoopArea) { TEST(degenerateLoop0) { // TODO: should we raise an error when given degenerate loops? - LatLng verts[] = {}; - _compareArea(verts, ARRAY_SIZE(verts), 0.0); + _compareArea(NULL, 0, 0.0); } } From 28105c3a51c833073bf22a16266d8abf14e088b8 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:20:59 -0500 Subject: [PATCH 48/74] remove docs --- website/docs/api/regions.mdx | 96 ------------------------------------ 1 file changed, 96 deletions(-) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 1431cc55bd..ceebf05ade 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -661,99 +661,3 @@ This function exists for memory management and is not exposed. - -## geoLoopAreaRads2 - -Compute the spherical area, in square radians, that is enclosed by a closed -`GeoLoop`. The loop must be a simple polygon traced counter-clockwise as viewed -from above the Earth, and each edge is interpreted as the shortest geodesic -between consecutive vertices (i.e., no side spans more than π radians). Areas -are normalized to the [0, 4π] interval, so reversing the vertex order yields -the complementary area of the loop. - - - - -```c -H3Error geoLoopAreaRads2(GeoLoop loop, double *out); -``` - -The caller supplies the `GeoLoop` memory (for example, by filling `loop.verts` -with previously allocated `LatLng` vertices) and receives the enclosed area in -square radians in `out`. Returns 0 (`E_SUCCESS`) on success. - - - - -:::note - -This function is not exposed in the Java bindings. - -::: - - - - -:::note - -This function is not exposed in the JavaScript bindings. - -::: - - - - -:::note - -This function isn't exposed directly, but can be used -via `h3.h3shape_area()`: - -```python -poly = h3.LatLngPoly([(90,0),(0,0),(0,90)]) -h3.h3shape_area(poly, unit='km^2') -``` - -::: - - - - -:::note - -This function is not exposed in the Go bindings. - -::: - - - - -:::note - -This function is not exposed in the DuckDB bindings. - -::: - - - - -:::note - -This function is not exposed in the CLI bindings. - -::: - - - - From a3ee94dac936611b64b2d7bbdb0a4579e2233234 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:22:57 -0500 Subject: [PATCH 49/74] drop H3_EXPORT --- src/apps/fuzzers/fuzzerGeoLoopArea.c | 2 +- src/apps/testapps/testGeoLoopArea.c | 4 ++-- src/h3lib/include/h3api.h.in | 2 +- src/h3lib/lib/area.c | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index b13b69f488..652c7865cf 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -28,7 +28,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { .verts = (LatLng *)data, }; double area; - H3_EXPORT(geoLoopAreaRads2)(loop, &area); + geoLoopAreaRads2(loop, &area); return 0; } diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index db331c71f8..dacdfa3265 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -31,7 +31,7 @@ static void _compareArea(LatLng *verts, int numVerts, double target_area) { GeoLoop loop = {.verts = verts, .numVerts = numVerts}; double out; - t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); + t_assertSuccess(geoLoopAreaRads2(loop, &out)); t_assert(fabs(out - target_area) < tol, "area should match"); } @@ -156,7 +156,7 @@ SUITE(geoLoopArea) { GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; double out; - t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); + t_assertSuccess(geoLoopAreaRads2(loop, &out)); double tol = 1e-13; if (t < 0.99) { diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index c8e7df2fee..bf13ff9d1c 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -355,7 +355,7 @@ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** * @brief Area in radians^2 enclosed by vertices of a GeoLoop. */ -DECLSPEC H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out); +DECLSPEC H3Error geoLoopAreaRads2(GeoLoop loop, double *out); /** @} */ /** @defgroup degsToRads degsToRads diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 5b0b2ac3cb..f1499357a5 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -76,7 +76,7 @@ static inline double cagnoli(LatLng x, LatLng y) { * @param out loop area in radians^2, in interval [0, 4*pi] * @return E_SUCCESS on success, or an error code otherwise */ -H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out) { +H3Error geoLoopAreaRads2(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms Adder adder = {0}; @@ -115,7 +115,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { } GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; - err = H3_EXPORT(geoLoopAreaRads2)(loop, out); + err = geoLoopAreaRads2(loop, out); if (NEVER(err)) { return err; } From 571d72ea4f9ac90e452ecd127c3d922b9d62c38d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:28:34 -0500 Subject: [PATCH 50/74] move geoLoopAreaRads2 out of public API for now --- CMakeLists.txt | 1 + src/apps/fuzzers/fuzzerGeoLoopArea.c | 1 + src/apps/testapps/testGeoLoopArea.c | 1 + src/h3lib/include/adder.h | 16 ++++++++++++++++ src/h3lib/include/area.h | 16 ++++++++++++++++ src/h3lib/include/h3api.h.in | 10 ---------- src/h3lib/lib/area.c | 2 ++ 7 files changed, 37 insertions(+), 10 deletions(-) create mode 100644 src/h3lib/include/area.h diff --git a/CMakeLists.txt b/CMakeLists.txt index d2bea788f9..2ee8421547 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -171,6 +171,7 @@ set(LIB_SOURCE_FILES src/h3lib/include/coordijk.h src/h3lib/include/algos.h src/h3lib/include/adder.h + src/h3lib/include/area.h src/h3lib/lib/h3Assert.c src/h3lib/lib/algos.c src/h3lib/lib/coordijk.c diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index 652c7865cf..eb39842b7a 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -18,6 +18,7 @@ */ #include "aflHarness.h" +#include "area.h" #include "h3api.h" #include "utility.h" diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index dacdfa3265..fd9ba1b230 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -22,6 +22,7 @@ #include +#include "area.h" #include "h3api.h" #include "test.h" #include "utility.h" diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 57682d4a78..dfebb0bcf1 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -1,3 +1,19 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + #ifndef ADDER_H #define ADDER_H diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h new file mode 100644 index 0000000000..22f43aeb2d --- /dev/null +++ b/src/h3lib/include/area.h @@ -0,0 +1,16 @@ +#ifndef AREA_H +#define AREA_H + +#include "h3api.h" + +/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 + * Functions for geoLoopAreaRads2 + * @{ + */ +/** + * @brief Area in radians^2 enclosed by vertices of a GeoLoop. + */ +H3Error geoLoopAreaRads2(GeoLoop loop, double *out); +/** @} */ + +#endif diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index bf13ff9d1c..f237dc3582 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,16 +348,6 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 - * Functions for geoLoopAreaRads2 - * @{ - */ -/** - * @brief Area in radians^2 enclosed by vertices of a GeoLoop. - */ -DECLSPEC H3Error geoLoopAreaRads2(GeoLoop loop, double *out); -/** @} */ - /** @defgroup degsToRads degsToRads * Functions for degsToRads * @{ diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index f1499357a5..e80d25da5f 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -18,6 +18,8 @@ * cells, polygons, multipolygons, etc.) */ +#include "area.h" + #include #include "adder.h" From 7d0f56a8488ccb447d6243bac80a925432f3173d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:31:23 -0500 Subject: [PATCH 51/74] slim down benchmark --- src/apps/benchmarks/benchmarkArea.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index fe1ffa4e29..193f6e66f2 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -27,7 +27,7 @@ static void doResSum(int res, bool print) { BEGIN_BENCHMARKS(); -static int MAX_RES = 4; +static int MAX_RES = 3; BENCHMARK(allCellsAtRes_print, 1, { for (int i = 0; i <= MAX_RES; i++) { @@ -35,7 +35,7 @@ BENCHMARK(allCellsAtRes_print, 1, { } }); -BENCHMARK(allCellsAtRes_noprint, 100, { +BENCHMARK(allCellsAtRes_noprint, 10, { for (int i = 0; i <= MAX_RES; i++) { doResSum(i, false); } From 77278a6089e1a83222c317fc78a87b0735ace5d6 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:50:26 -0500 Subject: [PATCH 52/74] clear these comments, maybe? --- src/h3lib/include/area.h | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index 22f43aeb2d..bf9235eb48 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -3,14 +3,6 @@ #include "h3api.h" -/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 - * Functions for geoLoopAreaRads2 - * @{ - */ -/** - * @brief Area in radians^2 enclosed by vertices of a GeoLoop. - */ H3Error geoLoopAreaRads2(GeoLoop loop, double *out); -/** @} */ #endif From fa555e33a2b063f84bcdf8347e266cd1db6ffabf Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:04:54 -0500 Subject: [PATCH 53/74] try adding area.h to APP_SOURCE_FILES --- CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ee8421547..f1a127a5ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,7 +200,8 @@ set(APP_SOURCE_FILES src/apps/applib/include/aflHarness.h src/apps/applib/lib/kml.c src/apps/applib/lib/utility.c - src/apps/applib/lib/args.c) + src/apps/applib/lib/args.c + src/h3lib/include/area.h) set(TEST_APP_SOURCE_FILES src/apps/applib/include/test.h src/apps/applib/lib/test.c) set(EXAMPLE_SOURCE_FILES From 64e9ad88115645e1301b7c74a0ade94cd9150742 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:14:12 -0500 Subject: [PATCH 54/74] can't fail if you don't try! Also, probably do not need to fuzz if no longer in the public API --- CMakeLists.txt | 5 +--- src/apps/fuzzers/README.md | 1 - src/apps/fuzzers/fuzzerGeoLoopArea.c | 37 ---------------------------- 3 files changed, 1 insertion(+), 42 deletions(-) delete mode 100644 src/apps/fuzzers/fuzzerGeoLoopArea.c diff --git a/CMakeLists.txt b/CMakeLists.txt index f1a127a5ab..6b4c8983ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,8 +200,7 @@ set(APP_SOURCE_FILES src/apps/applib/include/aflHarness.h src/apps/applib/lib/kml.c src/apps/applib/lib/utility.c - src/apps/applib/lib/args.c - src/h3lib/include/area.h) + src/apps/applib/lib/args.c) set(TEST_APP_SOURCE_FILES src/apps/applib/include/test.h src/apps/applib/lib/test.c) set(EXAMPLE_SOURCE_FILES @@ -312,7 +311,6 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerCellToChildPos.c src/apps/fuzzers/fuzzerInternalAlgos.c src/apps/fuzzers/fuzzerInternalCoordIjk.c - src/apps/fuzzers/fuzzerGeoLoopArea.c src/apps/benchmarks/benchmarkPolygonToCells.c src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c src/apps/benchmarks/benchmarkPolygon.c @@ -642,7 +640,6 @@ if(BUILD_FUZZERS) add_h3_fuzzer(fuzzerPolygonToCellsExperimentalNoHoles src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c) add_h3_fuzzer(fuzzerCellToChildPos src/apps/fuzzers/fuzzerCellToChildPos.c) - add_h3_fuzzer(fuzzerGeoLoopArea src/apps/fuzzers/fuzzerGeoLoopArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_fuzzer(fuzzerInternalAlgos src/apps/fuzzers/fuzzerInternalAlgos.c) diff --git a/src/apps/fuzzers/README.md b/src/apps/fuzzers/README.md index a559075950..192d4a37d5 100644 --- a/src/apps/fuzzers/README.md +++ b/src/apps/fuzzers/README.md @@ -65,7 +65,6 @@ The public API of H3 is covered in the following fuzzers: | originToDirectedEdges | [fuzzerDirectedEdge](./fuzzerDirectedEdge.c) | polygonToCells | [fuzzerPoylgonToCells](./fuzzerPolygonToCells.c) | polygonToCellsExperimental | [fuzzerPoylgonToCellsExperimental](./fuzzerPolygonToCellsExperimental.c) [fuzzerPoylgonToCellsExperimentalNoHoles](./fuzzerPolygonToCellsExperimentalNoHoles.c) -| geoLoopArea | [fuzzerGeoLoopArea](./fuzzerGeoLoopArea.c) | radsToDegs | Trivial | stringToH3 | [fuzzerIndexIO](./fuzzerIndexIO.c) | uncompactCells | [fuzzerCompact](./fuzzerCompact.c) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c deleted file mode 100644 index eb39842b7a..0000000000 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ /dev/null @@ -1,37 +0,0 @@ -/* - * Copyright 2025 Uber Technologies, Inc. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -/** @file - * @brief Fuzzer program for geoLoopAreaRads2 - */ - -#include "aflHarness.h" -#include "area.h" -#include "h3api.h" -#include "utility.h" - -int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { - // Interpret the input buffer as an array of LatLng vertices. - GeoLoop loop = { - .numVerts = (int)(size / sizeof(LatLng)), - .verts = (LatLng *)data, - }; - double area; - geoLoopAreaRads2(loop, &area); - - return 0; -} - -AFL_HARNESS_MAIN(sizeof(LatLng) * 1024); From a05939767c393eec9c04c0bfeb52b668a51819b0 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:28:15 -0500 Subject: [PATCH 55/74] minor --- CHANGELOG.md | 2 +- website/docs/api/regions.mdx | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e68810554..73f71e5dee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,7 @@ The public API of this library consists of the functions declared in file ## [Unreleased] ### Added - `reverseDirectedEdge` function (#1098) -- `geoLoopArea` function (#1101) +- (internal) `geoLoopArea` function (#1101) ### Changed - `cellAreaRads2` uses `geoLoopArea` (#1101) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index ceebf05ade..a5bf4c8510 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -661,3 +661,5 @@ This function exists for memory management and is not exposed. + + From b24442f334cd7552348025ed9a627304b67b3a78 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:54:42 -0500 Subject: [PATCH 56/74] remove justfile --- justfile | 53 ----------------------------------------------------- 1 file changed, 53 deletions(-) delete mode 100644 justfile diff --git a/justfile b/justfile deleted file mode 100644 index cffb5a6fde..0000000000 --- a/justfile +++ /dev/null @@ -1,53 +0,0 @@ -# TODO: Remove before landing - -init: purge - mkdir build - -build: - cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - -profile: init - cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make - # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - xcrun xctrace record \ - --template 'Time Profiler' \ - --output h3-prof.trace \ - --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon - open h3-prof.trace - -purge: - rm -rf build - rm -rf *.trace - rm -rf .ipynb_checkpoints - -test: build - # ./build/bin/testH3CellAreaExhaustive - # ./build/bin/testEdgeCellsToPoly - # ./build/bin/testDirectedEdge - # ./build/bin/testArea - # just test-slow - just test-fast - # ./build/bin/testGeoLoopArea - -time: - time ./build/bin/testArea - -test-fast: build - cd build; make test-fast - -test-slow: build - cd build; make test - -bench: build - ./build/bin/benchmarkArea - # ./build/bin/benchmarkCellsToLinkedMultiPolygon - # ./build/bin/benchmarkCellsToMultiPolygon - # ./build/bin/benchmarkDirectedEdge - -fuzz: - cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers - -coverage: purge - mkdir build - cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage - open build/coverage/index.html From dd316e0a1b30c647014118304293f42131641508 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 04:43:40 -0500 Subject: [PATCH 57/74] Adder adder = {}; --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/apps/testapps/testH3CellAreaExhaustive.c | 2 +- src/h3lib/include/adder.h | 2 +- src/h3lib/lib/area.c | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 193f6e66f2..cbeba2d8b8 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,7 +9,7 @@ #include "iterators.h" static void doResSum(int res, bool print) { - Adder adder = {0}; + Adder adder = {}; double cellArea; IterCellsResolution iter = iterInitRes(res); diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 5a8bd30aa6..784aa51bde 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -137,7 +137,7 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - Adder adder = {0}; + Adder adder = {}; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index dfebb0bcf1..9796ae520a 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -47,7 +47,7 @@ tradeoffs between Kahan, Neumaier, and `fsum`. ## Usage -Adder adder = {0}; // initialize +Adder adder = {}; // initialize kadd(&adder, x); // add x to summation out = adder.sum; // extract final sum */ diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index e80d25da5f..2394373e72 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -81,7 +81,7 @@ static inline double cagnoli(LatLng x, LatLng y) { H3Error geoLoopAreaRads2(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms - Adder adder = {0}; + Adder adder = {}; for (int i = 0; i < loop.numVerts; i++) { int j = (i + 1) % loop.numVerts; From c7275431a19cfa6b020774f612f7d205c40bd75a Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 04:41:17 -0500 Subject: [PATCH 58/74] (multi)polygon area functions --- justfile | 53 ++++++++++++++++++++++++++++++++++++++++ src/h3lib/include/area.h | 2 ++ src/h3lib/lib/area.c | 33 +++++++++++++++++++++++++ 3 files changed, 88 insertions(+) create mode 100644 justfile diff --git a/justfile b/justfile new file mode 100644 index 0000000000..cffb5a6fde --- /dev/null +++ b/justfile @@ -0,0 +1,53 @@ +# TODO: Remove before landing + +init: purge + mkdir build + +build: + cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + +profile: init + cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make + # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + xcrun xctrace record \ + --template 'Time Profiler' \ + --output h3-prof.trace \ + --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon + open h3-prof.trace + +purge: + rm -rf build + rm -rf *.trace + rm -rf .ipynb_checkpoints + +test: build + # ./build/bin/testH3CellAreaExhaustive + # ./build/bin/testEdgeCellsToPoly + # ./build/bin/testDirectedEdge + # ./build/bin/testArea + # just test-slow + just test-fast + # ./build/bin/testGeoLoopArea + +time: + time ./build/bin/testArea + +test-fast: build + cd build; make test-fast + +test-slow: build + cd build; make test + +bench: build + ./build/bin/benchmarkArea + # ./build/bin/benchmarkCellsToLinkedMultiPolygon + # ./build/bin/benchmarkCellsToMultiPolygon + # ./build/bin/benchmarkDirectedEdge + +fuzz: + cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers + +coverage: purge + mkdir build + cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage + open build/coverage/index.html diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index bf9235eb48..6d765bf440 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -4,5 +4,7 @@ #include "h3api.h" H3Error geoLoopAreaRads2(GeoLoop loop, double *out); +H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out); +H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out); #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 2394373e72..6b0e2e9f1d 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -124,3 +124,36 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { return E_SUCCESS; } + +H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out) { + Adder adder = {}; + double term; + + geoLoopAreaRads2(poly.geoloop, &term); + kadd(&adder, term); + + for (int i = 0; i < poly.numHoles; i++) { + geoLoopAreaRads2(poly.holes[i], &term); + + kadd(&adder, term); + kadd(&adder, -4.0 * M_PI); + } + + *out = adder.sum; + + return E_SUCCESS; +} + +H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) { + Adder adder = {}; + double term; + + for (int i = 0; i < mpoly.numPolygons; i++) { + geoPolygonAreaRads2(mpoly.polygons[i], &term); + kadd(&adder, term); + } + + *out = adder.sum; + + return E_SUCCESS; +} From 8c12e2f6fda1b9cf67730cfd8b7f4fd57716aa65 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:41:12 -0500 Subject: [PATCH 59/74] degenerate loop note --- src/apps/testapps/testGeoLoopArea.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index fd9ba1b230..a8331451d6 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -202,7 +202,8 @@ SUITE(geoLoopArea) { } TEST(degenerateLoop2) { - // TODO: should we raise an error when given degenerate loops? + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. LatLng verts[] = { {M_PI_2, 0.0}, {0.0, -M_PI_2}, @@ -211,7 +212,8 @@ SUITE(geoLoopArea) { } TEST(degenerateLoop1) { - // TODO: should we raise an error when given degenerate loops? + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. LatLng verts[] = { {0.0, 0.0}, }; @@ -219,7 +221,8 @@ SUITE(geoLoopArea) { } TEST(degenerateLoop0) { - // TODO: should we raise an error when given degenerate loops? + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. _compareArea(NULL, 0, 0.0); } } From 861bc8172d4884172d8339e026e58de7443beb69 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:03:56 -0500 Subject: [PATCH 60/74] destroy sigs --- src/h3lib/include/h3api.h.in | 9 +++++++++ src/h3lib/lib/algos.c | 4 ++++ 2 files changed, 13 insertions(+) diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index f237dc3582..de0e41d028 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -346,6 +346,15 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, /** @brief Free all memory created for a LinkedGeoPolygon */ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); + +/** @brief Free all memory created for a GeoLoop */ +DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop); + +/** @brief Free all memory created for a GeoPolygon */ +DECLSPEC void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly); + +/** @brief Free all memory created for a GeoMultiPolygon */ +DECLSPEC void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly); /** @} */ /** @defgroup degsToRads degsToRads diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index f21da3b3c6..f9b20c0aa5 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1280,3 +1280,7 @@ H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, } return normalizeResult; } + +void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) {} +void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) {} +void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly) {} From 59fefdcce47e7fea0e888e9c2db4c95dd052c635 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:14:41 -0500 Subject: [PATCH 61/74] destroy implementations --- src/h3lib/lib/algos.c | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index f9b20c0aa5..b46323f84f 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1281,6 +1281,31 @@ H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, return normalizeResult; } -void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) {} -void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) {} -void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly) {} +void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { + H3_MEMORY(free)(loop->verts); + loop->verts = NULL; + loop->numVerts = 0; +} + +void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) { + destroyGeoLoop(&poly->geoloop); + for (int i = 0; i < poly->numHoles; i++) { + destroyGeoLoop(&poly->holes[i]); + } + H3_MEMORY(free)(poly->holes); + poly->holes = NULL; + poly->numHoles = 0; +} + +/** + * Free all allocated memory for a GeoMultiPolygon. The caller is + * responsible for freeing memory allocated to input GeoMultiPolygon struct. + */ +void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly) { + for (int i = 0; i < mpoly->numPolygons; i++) { + destroyGeoPolygon(&mpoly->polygons[i]); + } + H3_MEMORY(free)(mpoly->polygons); + mpoly->polygons = NULL; + mpoly->numPolygons = 0; +} From 246f2f4fceafaf370652e1af9831a0c1dc4b6ff4 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:23:23 -0500 Subject: [PATCH 62/74] createGlobalMultiPolygon --- src/h3lib/include/area.h | 1 + src/h3lib/lib/area.c | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index 6d765bf440..962969ca5d 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -3,6 +3,7 @@ #include "h3api.h" +GeoMultiPolygon createGlobalMultiPolygon(); H3Error geoLoopAreaRads2(GeoLoop loop, double *out); H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out); H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out); diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 6b0e2e9f1d..87990c72d7 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -23,6 +23,7 @@ #include #include "adder.h" +#include "alloc.h" #include "constants.h" #include "h3Assert.h" #include "h3api.h" @@ -157,3 +158,37 @@ H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) { return E_SUCCESS; } + +GeoMultiPolygon createGlobalMultiPolygon() { + const int numPolygons = 8; + const int numVerts = 3; + const LatLng verts[numPolygons][numVerts] = { + {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, + {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, + {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, + {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, + {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, + {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, + {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, + {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, + }; + + GeoMultiPolygon mpoly = { + .numPolygons = numPolygons, + .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), + }; + + for (int i = 0; i < numPolygons; i++) { + GeoPolygon *poly = &mpoly.polygons[i]; + poly->numHoles = 0; + poly->holes = NULL; + poly->geoloop.numVerts = numVerts; + poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); + + for (int j = 0; j < numVerts; j++) { + poly->geoloop.verts[j] = verts[i][j]; + } + } + + return mpoly; +} From 4164d377c2d8c899a74bf36dc0786c33275e7147 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:35:20 -0500 Subject: [PATCH 63/74] testGeoMultiPolygon.c --- CMakeLists.txt | 1 + CMakeTests.cmake | 1 + justfile | 1 + src/apps/testapps/testGeoMultiPolygon.c | 41 +++++++++++++++++++++++++ 4 files changed, 44 insertions(+) create mode 100644 src/apps/testapps/testGeoMultiPolygon.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b4c8983ec..4f88c85551 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -282,6 +282,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c src/apps/testapps/testGeoLoopArea.c + src/apps/testapps/testGeoMultiPolygon.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 3508ff23a3..f7a2db4f26 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -260,6 +260,7 @@ add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) add_h3_test(testGeoLoopArea src/apps/testapps/testGeoLoopArea.c) +add_h3_test(testGeoMultiPolygon src/apps/testapps/testGeoMultiPolygon.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index cffb5a6fde..48362106d6 100644 --- a/justfile +++ b/justfile @@ -28,6 +28,7 @@ test: build # just test-slow just test-fast # ./build/bin/testGeoLoopArea + # ./build/bin/testGeoMultiPolygon time: time ./build/bin/testArea diff --git a/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c new file mode 100644 index 0000000000..9e42bb0466 --- /dev/null +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -0,0 +1,41 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** @file + * @brief tests for GeoMultiPolygon, GeoPolygon, and GeoLoop + * + * usage: `testGeoMultiPolygon` + */ + +#include + +#include "area.h" +#include "h3api.h" +#include "test.h" +#include "utility.h" + +SUITE(geoMultiPolygon) { + TEST(globalMultiPolygonArea) { + double tol = 1e-14; + double out; + + GeoMultiPolygon mpoly = createGlobalMultiPolygon(); + t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); + t_assert(fabs(out - 4 * M_PI) < tol, "area should match"); + + destroyGeoMultiPolygon(&mpoly); + } +} From 808c398d6ec9c7d92b4b29414add818194f624a4 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:25:53 -0600 Subject: [PATCH 64/74] try constants --- src/h3lib/lib/area.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 87990c72d7..7f5f6e05ef 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -162,7 +162,7 @@ H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) { GeoMultiPolygon createGlobalMultiPolygon() { const int numPolygons = 8; const int numVerts = 3; - const LatLng verts[numPolygons][numVerts] = { + const LatLng verts[8][3] = { {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, From 8f1c085c3c0e63dd2213974b69cde0f7169735eb Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:29:36 -0600 Subject: [PATCH 65/74] H3_EXPORT --- src/apps/testapps/testGeoMultiPolygon.c | 2 +- src/h3lib/lib/algos.c | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c index 9e42bb0466..cb6bb3e1ab 100644 --- a/src/apps/testapps/testGeoMultiPolygon.c +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -36,6 +36,6 @@ SUITE(geoMultiPolygon) { t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); t_assert(fabs(out - 4 * M_PI) < tol, "area should match"); - destroyGeoMultiPolygon(&mpoly); + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); } } diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index b46323f84f..7415a5ffdb 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1288,9 +1288,9 @@ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { } void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) { - destroyGeoLoop(&poly->geoloop); + H3_EXPORT(destroyGeoLoop)(&poly->geoloop); for (int i = 0; i < poly->numHoles; i++) { - destroyGeoLoop(&poly->holes[i]); + H3_EXPORT(destroyGeoLoop)(&poly->holes[i]); } H3_MEMORY(free)(poly->holes); poly->holes = NULL; @@ -1303,7 +1303,7 @@ void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) { */ void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly) { for (int i = 0; i < mpoly->numPolygons; i++) { - destroyGeoPolygon(&mpoly->polygons[i]); + H3_EXPORT(destroyGeoPolygon)(&mpoly->polygons[i]); } H3_MEMORY(free)(mpoly->polygons); mpoly->polygons = NULL; From e6bbf551aa63eed39eb04725b803e221b65bce24 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:56:45 -0600 Subject: [PATCH 66/74] rename to createGlobeMultiPolygon --- src/apps/testapps/testGeoMultiPolygon.c | 2 +- src/h3lib/include/area.h | 2 +- src/h3lib/lib/area.c | 10 +++++++++- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c index cb6bb3e1ab..60a153ffd2 100644 --- a/src/apps/testapps/testGeoMultiPolygon.c +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -32,7 +32,7 @@ SUITE(geoMultiPolygon) { double tol = 1e-14; double out; - GeoMultiPolygon mpoly = createGlobalMultiPolygon(); + GeoMultiPolygon mpoly = createGlobeMultiPolygon(); t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); t_assert(fabs(out - 4 * M_PI) < tol, "area should match"); diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index 962969ca5d..e8e17fa481 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -3,7 +3,7 @@ #include "h3api.h" -GeoMultiPolygon createGlobalMultiPolygon(); +GeoMultiPolygon createGlobeMultiPolygon(); H3Error geoLoopAreaRads2(GeoLoop loop, double *out); H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out); H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out); diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 7f5f6e05ef..fb9461c6ce 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -159,7 +159,15 @@ H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) { return E_SUCCESS; } -GeoMultiPolygon createGlobalMultiPolygon() { +/** + * Allocate a GeoMultiPolygon representing the entire globe. + * The globe is represented using 8 triangular polygons, with + * all edge arcs of exactly 90 degrees (i.e., pi/2 radians). + * Memory should be freed with `destroyGeoMultiPolygon`. + * + * @return GeoMultiPolygon covering entire globe + */ +GeoMultiPolygon createGlobeMultiPolygon() { const int numPolygons = 8; const int numVerts = 3; const LatLng verts[8][3] = { From c24fffa255a73921474f57d171dde54932cbabda Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 09:01:14 -0600 Subject: [PATCH 67/74] handle errors. add function docs --- src/h3lib/lib/algos.c | 8 ++++++++ src/h3lib/lib/area.c | 40 +++++++++++++++++++++++++++++++++++++--- 2 files changed, 45 insertions(+), 3 deletions(-) diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index 7415a5ffdb..ca8e54de9c 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1281,12 +1281,20 @@ H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, return normalizeResult; } +/** + * Free all allocated memory for a GeoLoop. The caller is + * responsible for freeing memory allocated to input GeoLoop struct. + */ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { H3_MEMORY(free)(loop->verts); loop->verts = NULL; loop->numVerts = 0; } +/** + * Free all allocated memory for a GeoPolygon. The caller is + * responsible for freeing memory allocated to input GeoPolygon struct. + */ void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) { H3_EXPORT(destroyGeoLoop)(&poly->geoloop); for (int i = 0; i < poly->numHoles; i++) { diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index fb9461c6ce..207bf7954b 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -126,16 +126,37 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { return E_SUCCESS; } +/** + * Area of GeoPolygon in radians^2. + * + * Outer GeoLoop vertices should be in counter-clockwise order. + * Hole GeoLoop vertices should be in clockwise order. + * Returned area is the area contained by the outer loop, minus + * the areas of the holes. See `geoLoopAreaRads2` for the expected + * form of GeoLoops. + * + * No check is made to ensure holes are disjoint or are contained + * within the outer GeoLoop. + * + * @param poly GeoPolygon + * @param out GeoPolygon area in radians^2 + * @return E_SUCCESS on success; error code otherwise + */ H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out) { + H3Error err; Adder adder = {}; double term; - geoLoopAreaRads2(poly.geoloop, &term); + err = geoLoopAreaRads2(poly.geoloop, &term); + if (err) return err; kadd(&adder, term); for (int i = 0; i < poly.numHoles; i++) { - geoLoopAreaRads2(poly.holes[i], &term); + err = geoLoopAreaRads2(poly.holes[i], &term); + if (err) return err; + // Due to clockwise order, holes will contribute area + // of "everything except the hole", so adjust with -4*pi term. kadd(&adder, term); kadd(&adder, -4.0 * M_PI); } @@ -145,12 +166,25 @@ H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out) { return E_SUCCESS; } +/** + * Area of GeoMultiPolygon in radians^2. + * + * Area is the sum of the areas of the polygons contained by `mpoly`. + * See `geoPolygonAreaRads2` for expected polygon format. + * No check is made to ensure polygons are disjoint. + * + * @param mpoly GeoMultiPolygon + * @param out GeoMultiPolygon area in radians^2 + * @return E_SUCCESS on success; error code otherwise + */ H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) { + H3Error err; Adder adder = {}; double term; for (int i = 0; i < mpoly.numPolygons; i++) { - geoPolygonAreaRads2(mpoly.polygons[i], &term); + err = geoPolygonAreaRads2(mpoly.polygons[i], &term); + if (err) return err; kadd(&adder, term); } From a154a82cbfa7af634fb6492c712e6ff81148372a Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 13:51:38 -0800 Subject: [PATCH 68/74] test for geoMultiPolygonAreaRads2 --- justfile | 4 +-- src/apps/testapps/testGeoMultiPolygon.c | 38 +++++++++++++++++++++++++ 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/justfile b/justfile index 48362106d6..23a76f6ed7 100644 --- a/justfile +++ b/justfile @@ -26,9 +26,9 @@ test: build # ./build/bin/testDirectedEdge # ./build/bin/testArea # just test-slow - just test-fast + # just test-fast # ./build/bin/testGeoLoopArea - # ./build/bin/testGeoMultiPolygon + ./build/bin/testGeoMultiPolygon time: time ./build/bin/testArea diff --git a/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c index 60a153ffd2..e033bc3a6c 100644 --- a/src/apps/testapps/testGeoMultiPolygon.c +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -38,4 +38,42 @@ SUITE(geoMultiPolygon) { H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); } + + TEST(triforceArea) { + LatLng _outer[] = { + {M_PI_2, 0}, + {0, 0}, + {0, M_PI_2}, + }; + LatLng _hole[] = { + {M_PI_2, 0}, + {0, M_PI_2}, + {0, 0}, + }; + + GeoLoop outer = { + .numVerts = 3, + .verts = _outer, + }; + GeoLoop hole = { + .numVerts = 3, + .verts = _hole, + }; + + GeoPolygon poly = { + .geoloop = outer, + .numHoles = 1, + .holes = &hole, + }; + + GeoMultiPolygon mpoly = { + .numPolygons = 1, + .polygons = &poly, + }; + + double tol = 1e-14; + double out; + t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); + t_assert(fabs(out - 0) < tol, "area should match"); + } } From fb79132d1578e9bd76ec4bed19fb4399b15bb772 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 14:27:20 -0800 Subject: [PATCH 69/74] better test; more coverage --- src/apps/testapps/testGeoMultiPolygon.c | 38 ++++++++++++++++++++----- 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c index e033bc3a6c..aedc22c126 100644 --- a/src/apps/testapps/testGeoMultiPolygon.c +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -21,7 +21,9 @@ */ #include +#include +#include "alloc.h" #include "area.h" #include "h3api.h" #include "test.h" @@ -39,13 +41,29 @@ SUITE(geoMultiPolygon) { H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); } - TEST(triforceArea) { + TEST(holeSameAsOuter) { + /** + * TODO: Replace with simpler test. + * + * I needed a test to exercize the "hole" branches of + * `destroyGeoMultiPolygon` and `geoMultiPolygonAreaRads2`. + * + * This is a verbose test because we have to allocate the + * `GeoMultiPolygon`. When we add in the `cellsToMultiPolygon` function, + * I can replace this with a much shorter, clearer test. + */ + + // Create a polygon with a triangle outer and a hole of exactly + // the same size, so the resulting polygon and multipolygons should + // have 0 area. LatLng _outer[] = { + // Counter-clockwise points {M_PI_2, 0}, {0, 0}, {0, M_PI_2}, }; LatLng _hole[] = { + // Same as above, but clockwise points {M_PI_2, 0}, {0, M_PI_2}, {0, 0}, @@ -53,27 +71,33 @@ SUITE(geoMultiPolygon) { GeoLoop outer = { .numVerts = 3, - .verts = _outer, + .verts = H3_MEMORY(malloc)(3 * sizeof(LatLng)), }; GeoLoop hole = { .numVerts = 3, - .verts = _hole, + .verts = H3_MEMORY(malloc)(3 * sizeof(LatLng)), }; + memcpy(outer.verts, _outer, 3 * sizeof(LatLng)); + memcpy(hole.verts, _hole, 3 * sizeof(LatLng)); + GeoPolygon poly = { .geoloop = outer, .numHoles = 1, - .holes = &hole, + .holes = H3_MEMORY(malloc)(sizeof(GeoLoop)), }; + poly.holes[0] = hole; GeoMultiPolygon mpoly = { .numPolygons = 1, - .polygons = &poly, + .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon)), }; + mpoly.polygons[0] = poly; - double tol = 1e-14; double out; t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); - t_assert(fabs(out - 0) < tol, "area should match"); + t_assert(fabs(out) < 1e-14, "Area should be 0"); + + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); } } From bedf3fae82674eb0ab23dd21d7e2349c4addf887 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 17 Dec 2025 15:33:47 -0800 Subject: [PATCH 70/74] fix missing brace from github merge --- justfile | 4 ++-- src/h3lib/lib/area.c | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/justfile b/justfile index 23a76f6ed7..47f436f8ed 100644 --- a/justfile +++ b/justfile @@ -25,10 +25,10 @@ test: build # ./build/bin/testEdgeCellsToPoly # ./build/bin/testDirectedEdge # ./build/bin/testArea - # just test-slow + just test-slow # just test-fast # ./build/bin/testGeoLoopArea - ./build/bin/testGeoMultiPolygon + # ./build/bin/testGeoMultiPolygon time: time ./build/bin/testArea diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 3f0d3caeeb..96aff2c844 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -233,6 +233,7 @@ GeoMultiPolygon createGlobeMultiPolygon() { } return mpoly; +} /** * Area of H3 cell in kilometers^2. From bdb05216d6377c6b043c3923d67f9908172d587a Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 17 Dec 2025 16:14:05 -0800 Subject: [PATCH 71/74] remove destroyGeoLoop and destroyGeoPolygon from public API --- src/h3lib/include/h3api.h.in | 6 ------ src/h3lib/lib/algos.c | 10 +++++----- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index de0e41d028..cdb9622f0b 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -347,12 +347,6 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, /** @brief Free all memory created for a LinkedGeoPolygon */ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); -/** @brief Free all memory created for a GeoLoop */ -DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop); - -/** @brief Free all memory created for a GeoPolygon */ -DECLSPEC void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly); - /** @brief Free all memory created for a GeoMultiPolygon */ DECLSPEC void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly); /** @} */ diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index ca8e54de9c..7a1d67e5e5 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1285,7 +1285,7 @@ H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, * Free all allocated memory for a GeoLoop. The caller is * responsible for freeing memory allocated to input GeoLoop struct. */ -void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { +void destroyGeoLoop(GeoLoop *loop) { H3_MEMORY(free)(loop->verts); loop->verts = NULL; loop->numVerts = 0; @@ -1295,10 +1295,10 @@ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { * Free all allocated memory for a GeoPolygon. The caller is * responsible for freeing memory allocated to input GeoPolygon struct. */ -void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) { - H3_EXPORT(destroyGeoLoop)(&poly->geoloop); +void destroyGeoPolygon(GeoPolygon *poly) { + destroyGeoLoop(&poly->geoloop); for (int i = 0; i < poly->numHoles; i++) { - H3_EXPORT(destroyGeoLoop)(&poly->holes[i]); + destroyGeoLoop(&poly->holes[i]); } H3_MEMORY(free)(poly->holes); poly->holes = NULL; @@ -1311,7 +1311,7 @@ void H3_EXPORT(destroyGeoPolygon)(GeoPolygon *poly) { */ void H3_EXPORT(destroyGeoMultiPolygon)(GeoMultiPolygon *mpoly) { for (int i = 0; i < mpoly->numPolygons; i++) { - H3_EXPORT(destroyGeoPolygon)(&mpoly->polygons[i]); + destroyGeoPolygon(&mpoly->polygons[i]); } H3_MEMORY(free)(mpoly->polygons); mpoly->polygons = NULL; From ff56f5a681f1d06b36b04cdefea530a1d5886928 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 17 Dec 2025 20:51:32 -0800 Subject: [PATCH 72/74] move global poly function --- justfile | 4 +-- src/apps/testapps/testGeoMultiPolygon.c | 1 + src/h3lib/include/algos.h | 3 ++ src/h3lib/include/area.h | 1 - src/h3lib/lib/algos.c | 42 +++++++++++++++++++++++++ src/h3lib/lib/area.c | 42 ------------------------- 6 files changed, 48 insertions(+), 45 deletions(-) diff --git a/justfile b/justfile index 47f436f8ed..48362106d6 100644 --- a/justfile +++ b/justfile @@ -25,8 +25,8 @@ test: build # ./build/bin/testEdgeCellsToPoly # ./build/bin/testDirectedEdge # ./build/bin/testArea - just test-slow - # just test-fast + # just test-slow + just test-fast # ./build/bin/testGeoLoopArea # ./build/bin/testGeoMultiPolygon diff --git a/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c index aedc22c126..5ea8bbcfba 100644 --- a/src/apps/testapps/testGeoMultiPolygon.c +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -23,6 +23,7 @@ #include #include +#include "algos.h" #include "alloc.h" #include "area.h" #include "h3api.h" diff --git a/src/h3lib/include/algos.h b/src/h3lib/include/algos.h index 0589c1ecc8..39808720f7 100644 --- a/src/h3lib/include/algos.h +++ b/src/h3lib/include/algos.h @@ -56,4 +56,7 @@ H3Error _gridDiskDistancesInternal(H3Index origin, int k, H3Index *out, // The safe gridRing algorithm. H3Error _gridRingInternal(H3Index origin, int k, H3Index *out); + +// Create a GeoMultiPolygon covering the entire globe +GeoMultiPolygon createGlobeMultiPolygon(); #endif diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index a4d3d8046b..1860f09bb9 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -22,7 +22,6 @@ #include "h3api.h" -GeoMultiPolygon createGlobeMultiPolygon(); H3Error geoLoopAreaRads2(GeoLoop loop, double *out); H3Error geoPolygonAreaRads2(GeoPolygon poly, double *out); H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out); diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c index 7a1d67e5e5..068b99e608 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1281,6 +1281,48 @@ H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, return normalizeResult; } +/** + * Allocate a GeoMultiPolygon representing the entire globe. + * The globe is represented using 8 triangular polygons, with + * all edge arcs of exactly 90 degrees (i.e., pi/2 radians). + * Memory should be freed with `destroyGeoMultiPolygon`. + * + * @return GeoMultiPolygon covering entire globe + */ +GeoMultiPolygon createGlobeMultiPolygon() { + const int numPolygons = 8; + const int numVerts = 3; + const LatLng verts[8][3] = { + {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, + {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, + {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, + {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, + {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, + {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, + {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, + {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, + }; + + GeoMultiPolygon mpoly = { + .numPolygons = numPolygons, + .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), + }; + + for (int i = 0; i < numPolygons; i++) { + GeoPolygon *poly = &mpoly.polygons[i]; + poly->numHoles = 0; + poly->holes = NULL; + poly->geoloop.numVerts = numVerts; + poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); + + for (int j = 0; j < numVerts; j++) { + poly->geoloop.verts[j] = verts[i][j]; + } + } + + return mpoly; +} + /** * Free all allocated memory for a GeoLoop. The caller is * responsible for freeing memory allocated to input GeoLoop struct. diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 96aff2c844..c8375cf48b 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -193,48 +193,6 @@ H3Error geoMultiPolygonAreaRads2(GeoMultiPolygon mpoly, double *out) { return E_SUCCESS; } -/** - * Allocate a GeoMultiPolygon representing the entire globe. - * The globe is represented using 8 triangular polygons, with - * all edge arcs of exactly 90 degrees (i.e., pi/2 radians). - * Memory should be freed with `destroyGeoMultiPolygon`. - * - * @return GeoMultiPolygon covering entire globe - */ -GeoMultiPolygon createGlobeMultiPolygon() { - const int numPolygons = 8; - const int numVerts = 3; - const LatLng verts[8][3] = { - {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, - {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, - {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, - {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, - {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, - {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, - {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, - {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, - }; - - GeoMultiPolygon mpoly = { - .numPolygons = numPolygons, - .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), - }; - - for (int i = 0; i < numPolygons; i++) { - GeoPolygon *poly = &mpoly.polygons[i]; - poly->numHoles = 0; - poly->holes = NULL; - poly->geoloop.numVerts = numVerts; - poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); - - for (int j = 0; j < numVerts; j++) { - poly->geoloop.verts[j] = verts[i][j]; - } - } - - return mpoly; -} - /** * Area of H3 cell in kilometers^2. * From 82fdb501c582a07491a8f5dbc263a55b1c995a0d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 17 Dec 2025 20:54:42 -0800 Subject: [PATCH 73/74] remove unnecessary import --- src/h3lib/lib/area.c | 1 - 1 file changed, 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index c8375cf48b..5836b063c9 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -23,7 +23,6 @@ #include #include "adder.h" -#include "alloc.h" #include "constants.h" #include "h3Assert.h" #include "h3api.h" From f65a3654146520f260084f1491260285bf980660 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 17 Dec 2025 21:11:04 -0800 Subject: [PATCH 74/74] remove justfile --- justfile | 54 ------------------------------------------------------ 1 file changed, 54 deletions(-) delete mode 100644 justfile diff --git a/justfile b/justfile deleted file mode 100644 index 48362106d6..0000000000 --- a/justfile +++ /dev/null @@ -1,54 +0,0 @@ -# TODO: Remove before landing - -init: purge - mkdir build - -build: - cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - -profile: init - cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make - # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - xcrun xctrace record \ - --template 'Time Profiler' \ - --output h3-prof.trace \ - --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon - open h3-prof.trace - -purge: - rm -rf build - rm -rf *.trace - rm -rf .ipynb_checkpoints - -test: build - # ./build/bin/testH3CellAreaExhaustive - # ./build/bin/testEdgeCellsToPoly - # ./build/bin/testDirectedEdge - # ./build/bin/testArea - # just test-slow - just test-fast - # ./build/bin/testGeoLoopArea - # ./build/bin/testGeoMultiPolygon - -time: - time ./build/bin/testArea - -test-fast: build - cd build; make test-fast - -test-slow: build - cd build; make test - -bench: build - ./build/bin/benchmarkArea - # ./build/bin/benchmarkCellsToLinkedMultiPolygon - # ./build/bin/benchmarkCellsToMultiPolygon - # ./build/bin/benchmarkDirectedEdge - -fuzz: - cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers - -coverage: purge - mkdir build - cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage - open build/coverage/index.html