diff --git a/src/apps/testapps/testGridPathCells.c b/src/apps/testapps/testGridPathCells.c index a99d33c19f..5d8b4c1834 100644 --- a/src/apps/testapps/testGridPathCells.c +++ b/src/apps/testapps/testGridPathCells.c @@ -19,6 +19,7 @@ * usage: `testGridPathCells` */ +#include #include #include #include @@ -29,6 +30,20 @@ #include "test.h" #include "utility.h" +static void assertPathValid(H3Index start, H3Index end, const H3Index *path, + int64_t size) { + t_assert(size > 0, "path size is positive"); + t_assert(path[0] == start, "path starts with start index"); + t_assert(path[size - 1] == end, "path ends with end index"); + + for (int64_t i = 1; i < size; i++) { + int isNeighbor; + t_assertSuccess( + H3_EXPORT(areNeighborCells)(path[i], path[i - 1], &isNeighbor)); + t_assert(isNeighbor, "path is contiguous"); + } +} + SUITE(gridPathCells) { TEST(gridPathCells_acrossMultipleFaces) { H3Index start = 0x85285aa7fffffff; @@ -40,15 +55,40 @@ SUITE(gridPathCells) { "Line not computable across multiple icosa faces"); } - TEST(gridPathCells_pentagon) { + TEST(gridPathCells_pentagonReverseInterpolation) { H3Index start = 0x820807fffffffff; H3Index end = 0x8208e7fffffffff; int64_t size; t_assertSuccess(H3_EXPORT(gridPathCellsSize)(start, end, &size)); - H3Index *path = calloc(sizeof(H3Index), size); - t_assert(H3_EXPORT(gridPathCells)(start, end, path) == E_PENTAGON, - "Line not computable due to pentagon distortion"); + H3Index *path = calloc(size, sizeof(H3Index)); + + t_assertSuccess(H3_EXPORT(gridPathCells)(start, end, path)); + assertPathValid(start, end, path, size); + + free(path); + } + + TEST(gridPathCells_knownFailureNotCoveredByReverseInterpolation) { + // Known limitation case: + // + // There are still pairs where `gridDistance` succeeds but interpolation + // fails in both origin-anchored local IJK charts (anchored at `start` + // and anchored at `end`). Since `gridPathCells` only performs these two + // interpolation attempts, it returns an error. + // + // This test keeps a pinned input pair to demonstrate that the current + // approach is not complete. + H3Index start = 0x8411b61ffffffff; + H3Index end = 0x84016d3ffffffff; + + int64_t size; + t_assertSuccess(H3_EXPORT(gridPathCellsSize)(start, end, &size)); + H3Index *path = calloc(size, sizeof(H3Index)); + + H3Error err = H3_EXPORT(gridPathCells)(start, end, path); + t_assert(err != E_SUCCESS, "Expected gridPathCells to fail"); + free(path); } } diff --git a/src/h3lib/lib/localij.c b/src/h3lib/lib/localij.c index 0f96adfa19..489023f034 100644 --- a/src/h3lib/lib/localij.c +++ b/src/h3lib/lib/localij.c @@ -22,9 +22,12 @@ #include #include #include +#include #include #include +#include "algos.h" +#include "alloc.h" #include "baseCells.h" #include "faceijk.h" #include "h3Assert.h" @@ -398,9 +401,10 @@ H3Error localIjkToCell(H3Index origin, const CoordIJK *ijk, H3Index *out) { for (int i = 0; i < pentagonRotations; i++) { dir = _rotate60ccw(dir); } - // The pentagon rotations are being chosen so that dir is not the - // deleted direction. If it still happens, it means we're moving - // into a deleted subsequence, so there is no index here. + // The pentagon rotations are chosen to avoid the deleted direction + // (the missing neighbor direction around a pentagon). If we still + // land on it, the coordinate would cross pentagon distortion and + // cannot be represented. if (dir == K_AXES_DIGIT) { return E_PENTAGON; } @@ -603,13 +607,15 @@ H3Error H3_EXPORT(gridDistance)(H3Index origin, H3Index index, int64_t *out) { /** * Number of indexes in a line from the start index to the end index, - * to be used for allocating memory. Returns a negative number if the - * line cannot be computed. + * to be used for allocating memory. + * + * On success, sets `*size` to `gridDistance(start, end) + 1` + * (including both endpoints). * * @param start Start index of the line * @param end End index of the line - * @param size Size of the line - * @returns 0 on success, or another value on error + * @param size Output size of the line + * @returns 0 on success, otherwise the error from `gridDistance` */ H3Error H3_EXPORT(gridPathCellsSize)(H3Index start, H3Index end, int64_t *size) { @@ -654,34 +660,32 @@ static void cubeRound(double i, double j, double k, CoordIJK *ijk) { } /** - * Given two H3 indexes, return the line of indexes between them (inclusive). - * - * This function may fail to find the line between two indexes, for - * example if they are very far apart. It may also fail when finding - * distances for indexes on opposite sides of a pentagon. + * Attempts to generate a shortest-length path by interpolating through an + * origin-anchored local IJK coordinate space. * - * Notes: + * This helper implements the interpolation-based path construction used by + * `gridPathCells`. It can fail if interpolation lands on intermediate IJK + * coordinates that cannot be mapped back to valid H3 cells. This can occur + * because the origin-anchored local IJ(K) space is not globally continuous, and + * some intermediate coordinates do not have an inverse mapping back to a cell + * in the chosen chart (for example, due to discontinuities or warping near + * pentagons). * - * - The specific output of this function should not be considered stable - * across library versions. The only guarantees the library provides are - * that the line length will be `gridDistance(start, end) + 1` and that - * every index in the line will be a neighbor of the preceding index. - * - Lines are drawn in grid space, and may not correspond exactly to either - * Cartesian lines or great arcs. + * The output is written to `out[outOffset + outStep * n]`, allowing callers to + * fill the path in either direction without an intermediate buffer. * - * @param start Start index of the line - * @param end End index of the line - * @param out Output array, which must be of size gridPathCellsSize(start, end) - * @return 0 on success, or another value on failure. + * @param start Origin cell for local IJK conversions + * @param end Target cell + * @param distance Expected edge distance between `start` and `end` + * @param out Output buffer + * @param outOffset Output index for the first element + * @param outStep Output stride (+1 for forward fill, -1 for reverse fill) + * @return E_SUCCESS if all intermediate steps convert successfully, otherwise + * the first encountered conversion error */ -H3Error H3_EXPORT(gridPathCells)(H3Index start, H3Index end, H3Index *out) { - int64_t distance; - H3Error distanceError = H3_EXPORT(gridDistance)(start, end, &distance); - // Early exit if we can't calculate the line - if (distanceError) { - return distanceError; - } - +static H3Error gridPathCellsInterpolate(H3Index start, H3Index end, + int64_t distance, H3Index *out, + int64_t outOffset, int64_t outStep) { // Get IJK coords for the start and end. We've already confirmed // that these can be calculated with the distance check above. CoordIJK startIjk = {0}; @@ -703,7 +707,7 @@ H3Error H3_EXPORT(gridPathCells)(H3Index start, H3Index end, H3Index *out) { ijkToCube(&startIjk); ijkToCube(&endIjk); - double invDistance = distance ? 1.0 / (double)distance : 0; + double invDistance = 1.0 / (double)distance; double iStep = (double)(endIjk.i - startIjk.i) * invDistance; double jStep = (double)(endIjk.j - startIjk.j) * invDistance; double kStep = (double)(endIjk.k - startIjk.k) * invDistance; @@ -715,7 +719,8 @@ H3Error H3_EXPORT(gridPathCells)(H3Index start, H3Index end, H3Index *out) { (double)startIjk.k + kStep * n, ¤tIjk); // Convert cube -> ijk -> h3 index cubeToIjk(¤tIjk); - H3Error currentError = localIjkToCell(start, ¤tIjk, &out[n]); + const int64_t idx = outOffset + outStep * n; + H3Error currentError = localIjkToCell(start, ¤tIjk, &out[idx]); if (currentError) { // The cells between `start` and `end` may fall in pentagon // distortion. @@ -725,3 +730,65 @@ H3Error H3_EXPORT(gridPathCells)(H3Index start, H3Index end, H3Index *out) { return E_SUCCESS; } + +/** + * Given two H3 indexes, return the line of indexes between them (inclusive). + * + * This function relies on `gridDistance(start, end)` to determine the expected + * path length, and will return the same error if `gridDistance` fails. + * + * Path construction is performed by straight-line interpolation in the + * origin-anchored local IJK coordinate space: + * + * - First, interpolate using `start` as the local IJK origin. + * - If interpolation fails, retry using `end` as the local IJK origin and + * reverse the resulting sequence into `out`. + * + * If both interpolation attempts fail, this function returns the error from the + * first attempt. + * + * Notes: + * + * - The specific output of this function should not be considered stable + * across library versions. The only guarantees the library provides are + * that the line length will be `gridDistance(start, end) + 1` and that + * every index in the line will be a neighbor of the preceding index. + * - Lines are drawn in grid space, and may not correspond exactly to either + * Cartesian lines or great arcs. + * + * @param start Start index of the line + * @param end End index of the line + * @param out Output array, which must be of size gridPathCellsSize(start, end) + * @return 0 on success, or another value on failure. + */ +H3Error H3_EXPORT(gridPathCells)(H3Index start, H3Index end, H3Index *out) { + int64_t distance; + H3Error distanceError = H3_EXPORT(gridDistance)(start, end, &distance); + // Early exit if we can't calculate the line + if (distanceError) { + return distanceError; + } + + if (distance == 0) { + out[0] = start; + return E_SUCCESS; + } + + // Straight-line interpolation in local IJK space anchored at `start`. + H3Error interpolateErr = + gridPathCellsInterpolate(start, end, distance, out, 0, 1); + if (!interpolateErr) { + return E_SUCCESS; + } + + // Retry interpolation anchored at `end` and reverse the output. + // This can resolve cases where the local IJK chart is discontinuous + // relative to one origin but not the other. + H3Error reverseErr = + gridPathCellsInterpolate(end, start, distance, out, distance, -1); + if (!reverseErr) { + return E_SUCCESS; + } + + return interpolateErr; +}