(Multi)polygon helper functions

CHANGELOG.md

@@ -7,7 +7,11 @@ The public API of this library consists of the functions declared in file
 
 ## [Unreleased]
 ### Added
-- `reverseDirectedEdge` function
+- `reverseDirectedEdge` function (#1098)
+- (internal) `geoLoopArea` function (#1101)
+
+### Changed
+- `cellAreaRads2` uses `geoLoopArea` (#1101)
 
 ## [4.4.1] - 2025-11-11
 ### Fixed

CMakeLists.txt

@@ -170,6 +170,8 @@ 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/include/area.h
     src/h3lib/lib/h3Assert.c
     src/h3lib/lib/algos.c
     src/h3lib/lib/coordijk.c
@@ -188,7 +190,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 +281,8 @@ 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/testapps/testGeoMultiPolygon.c
     src/apps/miscapps/cellToBoundaryHier.c
     src/apps/miscapps/cellToLatLngHier.c
     src/apps/miscapps/generateBaseCellNeighbors.c
@@ -317,7 +322,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}
@@ -672,6 +678,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)

CMakeTests.cmake

@@ -259,6 +259,8 @@ 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(testGeoMultiPolygon src/apps/testapps/testGeoMultiPolygon.c)
 
 add_h3_test_with_arg(testH3NeighborRotations
                      src/apps/testapps/testH3NeighborRotations.c 0)

justfile

@@ -0,0 +1,54 @@
+# 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

src/apps/benchmarks/benchmarkArea.c

@@ -0,0 +1,44 @@
+#include <math.h>
+#include <stdbool.h>
+#include <stdio.h>
+
+#include "adder.h"
+#include "benchmark.h"
+#include "constants.h"
+#include "h3api.h"
+#include "iterators.h"
+
+static void doResSum(int res, bool print) {
+    Adder adder = {};
+    double cellArea;
+    IterCellsResolution iter = iterInitRes(res);
+
+    for (; iter.h; iterStepRes(&iter)) {
+        H3_EXPORT(cellAreaRads2)(iter.h, &cellArea);
+        kadd(&adder, cellArea);
+    }
+
+    double total_area = adder.sum;
+    double diff = fabs(total_area - 4 * M_PI);
+    if (print) {
+        printf("res: %d, diff: %e\n", res, diff);
+    }
+}
+
+BEGIN_BENCHMARKS();
+
+static int MAX_RES = 3;
+
+BENCHMARK(allCellsAtRes_print, 1, {
+    for (int i = 0; i <= MAX_RES; i++) {
+        doResSum(i, true);
+    }
+});
+
+BENCHMARK(allCellsAtRes_noprint, 10, {
+    for (int i = 0; i <= MAX_RES; i++) {
+        doResSum(i, false);
+    }
+});
+
+END_BENCHMARKS();

src/apps/testapps/testGeoLoopArea.c

@@ -0,0 +1,228 @@
+/*
+ * 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 H3 GeoLoop area calculation
+ *
+ *  usage: `testGeoLoopArea`
+ */
+
+#include <math.h>
+
+#include "area.h"
+#include "h3api.h"
+#include "test.h"
+#include "utility.h"
+
+static void _compareArea(LatLng *verts, int numVerts, double target_area) {
+    double tol = 1e-14;
+    GeoLoop loop = {.verts = verts, .numVerts = numVerts};
+
+    double out;
+    t_assertSuccess(geoLoopAreaRads2(loop, &out));
+    t_assert(fabs(out - target_area) < tol, "area should match");
+}
+
+SUITE(geoLoopArea) {
+    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},
+            {0.0, M_PI_2},
+        };
+
+        _compareArea(verts, ARRAY_SIZE(verts), M_PI / 2);
+    }
+
+    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},
+            {M_PI_2, 0.0},
+        };
+
+        _compareArea(verts, ARRAY_SIZE(verts), 7 * M_PI / 2);
+    }
+
+    TEST(slice) {
+        /*
+        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},
+            {-M_PI_2, 0.0},
+            {0.0, M_PI_2},
+        };
+
+        _compareArea(verts, ARRAY_SIZE(verts), M_PI);
+    }
+
+    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},
+            {-M_PI_2, 0.0},
+            {0.0, 0.0},
+        };
+
+        _compareArea(verts, ARRAY_SIZE(verts), 3 * M_PI);
+    }
+
+    TEST(hemisphere_east) {
+        /*
+        Stitch 4 1/8 triangles together to cover the eastern hemisphere.
+        */
+        LatLng verts[] = {
+            {M_PI_2, 0.0},
+            {0.0, 0},
+            {-M_PI_2, 0.0},
+            {0.0, M_PI},
+        };
+
+        _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI);
+    }
+
+    TEST(hemisphere_north) {
+        /*
+        Stitch 4 1/8 triangles together to cover the northern hemisphere.
+        */
+        LatLng verts[] = {
+            {0.0, -M_PI},
+            {0.0, -M_PI_2},
+            {0.0, 0.0},
+            {0.0, M_PI_2},
+        };
+
+        _compareArea(verts, ARRAY_SIZE(verts), 2 * 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},
+                {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(geoLoopAreaRads2(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).
+            */
+        }
+    }
+
+    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);
+    }
+
+    TEST(degenerateLoop2) {
+        // Note that `geoLoopAreaRads2()` works without error on degenerate
+        // loops, returning 0 area.
+        LatLng verts[] = {
+            {M_PI_2, 0.0},
+            {0.0, -M_PI_2},
+        };
+        _compareArea(verts, ARRAY_SIZE(verts), 0.0);
+    }
+
+    TEST(degenerateLoop1) {
+        // Note that `geoLoopAreaRads2()` works without error on degenerate
+        // loops, returning 0 area.
+        LatLng verts[] = {
+            {0.0, 0.0},
+        };
+        _compareArea(verts, ARRAY_SIZE(verts), 0.0);
+    }
+
+    TEST(degenerateLoop0) {
+        // Note that `geoLoopAreaRads2()` works without error on degenerate
+        // loops, returning 0 area.
+        _compareArea(NULL, 0, 0.0);
+    }
+}