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/src/apps/testapps/testGeoMultiPolygon.c b/src/apps/testapps/testGeoMultiPolygon.c new file mode 100644 index 0000000000..5ea8bbcfba --- /dev/null +++ b/src/apps/testapps/testGeoMultiPolygon.c @@ -0,0 +1,104 @@ +/* + * 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 + +#include "algos.h" +#include "alloc.h" +#include "area.h" +#include "h3api.h" +#include "test.h" +#include "utility.h" + +SUITE(geoMultiPolygon) { + TEST(globalMultiPolygonArea) { + double tol = 1e-14; + double out; + + GeoMultiPolygon mpoly = createGlobeMultiPolygon(); + t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); + t_assert(fabs(out - 4 * M_PI) < tol, "area should match"); + + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + } + + 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}, + }; + + GeoLoop outer = { + .numVerts = 3, + .verts = H3_MEMORY(malloc)(3 * sizeof(LatLng)), + }; + GeoLoop hole = { + .numVerts = 3, + .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 = H3_MEMORY(malloc)(sizeof(GeoLoop)), + }; + poly.holes[0] = hole; + + GeoMultiPolygon mpoly = { + .numPolygons = 1, + .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon)), + }; + mpoly.polygons[0] = poly; + + double out; + t_assertSuccess(geoMultiPolygonAreaRads2(mpoly, &out)); + t_assert(fabs(out) < 1e-14, "Area should be 0"); + + H3_EXPORT(destroyGeoMultiPolygon)(&mpoly); + } +} 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 c561f53773..1860f09bb9 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -23,5 +23,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/include/h3api.h.in b/src/h3lib/include/h3api.h.in index f237dc3582..cdb9622f0b 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -346,6 +346,9 @@ 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 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..068b99e608 100644 --- a/src/h3lib/lib/algos.c +++ b/src/h3lib/lib/algos.c @@ -1280,3 +1280,82 @@ 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. + */ +void 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 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; +} diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 166059fcd9..5836b063c9 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -125,6 +125,73 @@ 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; + + err = geoLoopAreaRads2(poly.geoloop, &term); + if (err) return err; + kadd(&adder, term); + + for (int i = 0; i < poly.numHoles; i++) { + 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); + } + + *out = adder.sum; + + 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++) { + err = geoPolygonAreaRads2(mpoly.polygons[i], &term); + if (err) return err; + kadd(&adder, term); + } + + *out = adder.sum; + + return E_SUCCESS; +} + /** * Area of H3 cell in kilometers^2. *