Coverage Report

Created: 2025-08-26 06:12

/src/h3/src/h3lib/lib/algos.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2016-2021 Uber Technologies, Inc.
3
 *
4
 * Licensed under the Apache License, Version 2.0 (the "License");
5
 * you may not use this file except in compliance with the License.
6
 * You may obtain a copy of the License at
7
 *
8
 *         http://www.apache.org/licenses/LICENSE-2.0
9
 *
10
 * Unless required by applicable law or agreed to in writing, software
11
 * distributed under the License is distributed on an "AS IS" BASIS,
12
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13
 * See the License for the specific language governing permissions and
14
 * limitations under the License.
15
 */
16
/** @file algos.c
17
 * @brief   Hexagon grid algorithms
18
 */
19
20
#include "algos.h"
21
22
#include <assert.h>
23
#include <float.h>
24
#include <math.h>
25
#include <stdbool.h>
26
#include <stdlib.h>
27
#include <string.h>
28
29
#include "alloc.h"
30
#include "baseCells.h"
31
#include "bbox.h"
32
#include "faceijk.h"
33
#include "h3Assert.h"
34
#include "h3Index.h"
35
#include "h3api.h"
36
#include "latLng.h"
37
#include "linkedGeo.h"
38
#include "polygon.h"
39
#include "vertexGraph.h"
40
41
/*
42
 * Return codes from gridDiskUnsafe and related functions.
43
 */
44
45
184M
#define MAX_ONE_RING_SIZE 7
46
9.35k
#define POLYGON_TO_CELLS_BUFFER 12
47
48
/**
49
 * Directions used for traversing a hexagonal ring counterclockwise around
50
 * {1, 0, 0}
51
 *
52
 * <pre>
53
 *      _
54
 *    _/ \\_
55
 *   / \\5/ \\
56
 *   \\0/ \\4/
57
 *   / \\_/ \\
58
 *   \\1/ \\3/
59
 *     \\2/
60
 * </pre>
61
 */
62
static const Direction DIRECTIONS[6] = {J_AXES_DIGIT, JK_AXES_DIGIT,
63
                                        K_AXES_DIGIT, IK_AXES_DIGIT,
64
                                        I_AXES_DIGIT, IJ_AXES_DIGIT};
65
66
/**
67
 * Direction used for traversing to the next outward hexagonal ring.
68
 */
69
static const Direction NEXT_RING_DIRECTION = I_AXES_DIGIT;
70
71
/**
72
 * New digit when traversing along class II grids.
73
 *
74
 * Current digit -> direction -> new digit.
75
 */
76
static const Direction NEW_DIGIT_II[7][7] = {
77
    {CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT,
78
     IK_AXES_DIGIT, IJ_AXES_DIGIT},
79
    {K_AXES_DIGIT, I_AXES_DIGIT, JK_AXES_DIGIT, IJ_AXES_DIGIT, IK_AXES_DIGIT,
80
     J_AXES_DIGIT, CENTER_DIGIT},
81
    {J_AXES_DIGIT, JK_AXES_DIGIT, K_AXES_DIGIT, I_AXES_DIGIT, IJ_AXES_DIGIT,
82
     CENTER_DIGIT, IK_AXES_DIGIT},
83
    {JK_AXES_DIGIT, IJ_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT, CENTER_DIGIT,
84
     K_AXES_DIGIT, J_AXES_DIGIT},
85
    {I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, J_AXES_DIGIT,
86
     JK_AXES_DIGIT, K_AXES_DIGIT},
87
    {IK_AXES_DIGIT, J_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, JK_AXES_DIGIT,
88
     IJ_AXES_DIGIT, I_AXES_DIGIT},
89
    {IJ_AXES_DIGIT, CENTER_DIGIT, IK_AXES_DIGIT, J_AXES_DIGIT, K_AXES_DIGIT,
90
     I_AXES_DIGIT, JK_AXES_DIGIT}};
91
92
/**
93
 * New traversal direction when traversing along class II grids.
94
 *
95
 * Current digit -> direction -> new ap7 move (at coarser level).
96
 */
97
static const Direction NEW_ADJUSTMENT_II[7][7] = {
98
    {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT,
99
     CENTER_DIGIT, CENTER_DIGIT},
100
    {CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT,
101
     IK_AXES_DIGIT, CENTER_DIGIT},
102
    {CENTER_DIGIT, CENTER_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
103
     CENTER_DIGIT, J_AXES_DIGIT},
104
    {CENTER_DIGIT, K_AXES_DIGIT, JK_AXES_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
105
     CENTER_DIGIT, CENTER_DIGIT},
106
    {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
107
     I_AXES_DIGIT, IJ_AXES_DIGIT},
108
    {CENTER_DIGIT, IK_AXES_DIGIT, CENTER_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
109
     IK_AXES_DIGIT, CENTER_DIGIT},
110
    {CENTER_DIGIT, CENTER_DIGIT, J_AXES_DIGIT, CENTER_DIGIT, IJ_AXES_DIGIT,
111
     CENTER_DIGIT, IJ_AXES_DIGIT}};
112
113
/**
114
 * New traversal direction when traversing along class III grids.
115
 *
116
 * Current digit -> direction -> new ap7 move (at coarser level).
117
 */
118
static const Direction NEW_DIGIT_III[7][7] = {
119
    {CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT,
120
     IK_AXES_DIGIT, IJ_AXES_DIGIT},
121
    {K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT,
122
     IJ_AXES_DIGIT, CENTER_DIGIT},
123
    {J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT,
124
     CENTER_DIGIT, K_AXES_DIGIT},
125
    {JK_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT,
126
     K_AXES_DIGIT, J_AXES_DIGIT},
127
    {I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT,
128
     J_AXES_DIGIT, JK_AXES_DIGIT},
129
    {IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT,
130
     JK_AXES_DIGIT, I_AXES_DIGIT},
131
    {IJ_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT,
132
     I_AXES_DIGIT, IK_AXES_DIGIT}};
133
134
/**
135
 * New traversal direction when traversing along class III grids.
136
 *
137
 * Current digit -> direction -> new ap7 move (at coarser level).
138
 */
139
static const Direction NEW_ADJUSTMENT_III[7][7] = {
140
    {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT,
141
     CENTER_DIGIT, CENTER_DIGIT},
142
    {CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
143
     K_AXES_DIGIT, CENTER_DIGIT},
144
    {CENTER_DIGIT, CENTER_DIGIT, J_AXES_DIGIT, J_AXES_DIGIT, CENTER_DIGIT,
145
     CENTER_DIGIT, IJ_AXES_DIGIT},
146
    {CENTER_DIGIT, JK_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
147
     CENTER_DIGIT, CENTER_DIGIT},
148
    {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
149
     IK_AXES_DIGIT, I_AXES_DIGIT},
150
    {CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT, CENTER_DIGIT, IK_AXES_DIGIT,
151
     IK_AXES_DIGIT, CENTER_DIGIT},
152
    {CENTER_DIGIT, CENTER_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
153
     CENTER_DIGIT, IJ_AXES_DIGIT}};
154
155
/**
156
 * k value which will encompass all cells at resolution 15.
157
 * This is the largest possible k in the H3 grid system.
158
 */
159
static const int K_ALL_CELLS_AT_RES_15 = 13780510;
160
161
/**
162
 * Maximum number of cells that result from the gridDisk algorithm with the
163
 * given k. Formula source and proof: https://oeis.org/A003215
164
 *
165
 * @param   k   k value, k >= 0.
166
 * @param out   size in indexes
167
 */
168
41.6k
H3Error H3_EXPORT(maxGridDiskSize)(int k, int64_t *out) {
169
41.6k
    if (k < 0) {
170
35
        return E_DOMAIN;
171
35
    }
172
41.6k
    if (k >= K_ALL_CELLS_AT_RES_15) {
173
        // If a k value of this value or above is provided, this function will
174
        // estimate more cells than exist in the H3 grid at the finest
175
        // resolution. This is a problem since the function does signed integer
176
        // arithmetic on `k`, which could overflow. To prevent that, instead
177
        // substitute the maximum number of cells in the grid, as it should not
178
        // be possible for the gridDisk functions to exceed that. Note this is
179
        // not resolution specific. So, when resolution < 15, this function may
180
        // still estimate a size larger than the number of cells in the grid.
181
21
        return H3_EXPORT(getNumCells)(MAX_H3_RES, out);
182
21
    }
183
41.6k
    *out = 3 * (int64_t)k * ((int64_t)k + 1) + 1;
184
41.6k
    return E_SUCCESS;
185
41.6k
}
186
187
/**
188
 * Produce cells within grid distance k of the origin cell.
189
 *
190
 * k-ring 0 is defined as the origin cell, k-ring 1 is defined as k-ring 0 and
191
 * all neighboring cells, and so on.
192
 *
193
 * Output is placed in the provided array in no particular order. Elements of
194
 * the output array may be left zero, as can happen when crossing a pentagon.
195
 *
196
 * @param  origin   origin cell
197
 * @param  k        k >= 0
198
 * @param  out      zero-filled array which must be of size maxGridDiskSize(k)
199
 */
200
23.1M
H3Error H3_EXPORT(gridDisk)(H3Index origin, int k, H3Index *out) {
201
23.1M
    return H3_EXPORT(gridDiskDistances)(origin, k, out, NULL);
202
23.1M
}
203
204
/**
205
 * Produce cells and their distances from the given origin cell, up to
206
 * distance k.
207
 *
208
 * k-ring 0 is defined as the origin cell, k-ring 1 is defined as k-ring 0 and
209
 * all neighboring cells, and so on.
210
 *
211
 * Output is placed in the provided array in no particular order. Elements of
212
 * the output array may be left zero, as can happen when crossing a pentagon.
213
 *
214
 * @param  origin      origin cell
215
 * @param  k           k >= 0
216
 * @param  out         zero-filled array which must be of size
217
 * maxGridDiskSize(k)
218
 * @param  distances   NULL or a zero-filled array which must be of size
219
 *                     maxGridDiskSize(k)
220
 */
221
H3Error H3_EXPORT(gridDiskDistances)(H3Index origin, int k, H3Index *out,
222
23.1M
                                     int *distances) {
223
    // Optimistically try the faster gridDiskUnsafe algorithm first
224
23.1M
    const H3Error failed =
225
23.1M
        H3_EXPORT(gridDiskDistancesUnsafe)(origin, k, out, distances);
226
23.1M
    if (failed) {
227
40.2k
        int64_t maxIdx;
228
40.2k
        H3Error err = H3_EXPORT(maxGridDiskSize)(k, &maxIdx);
229
40.2k
        if (err) {
230
0
            return err;
231
0
        }
232
        // Fast algo failed, fall back to slower, correct algo
233
        // and also wipe out array because contents untrustworthy
234
40.2k
        memset(out, 0, maxIdx * sizeof(H3Index));
235
236
40.2k
        if (distances == NULL) {
237
40.0k
            distances = H3_MEMORY(calloc)(maxIdx, sizeof(int));
238
40.0k
            if (!distances) {
239
0
                return E_MEMORY_ALLOC;
240
0
            }
241
40.0k
            H3Error result = _gridDiskDistancesInternal(origin, k, out,
242
40.0k
                                                        distances, maxIdx, 0);
243
40.0k
            H3_MEMORY(free)(distances);
244
40.0k
            return result;
245
40.0k
        } else {
246
242
            memset(distances, 0, maxIdx * sizeof(int));
247
242
            return _gridDiskDistancesInternal(origin, k, out, distances, maxIdx,
248
242
                                              0);
249
242
        }
250
23.0M
    } else {
251
23.0M
        return E_SUCCESS;
252
23.0M
    }
253
23.1M
}
254
255
/**
256
 * Internal algorithm for the safe but slow version of gridDiskDistances
257
 *
258
 * Adds the origin cell to the output set (treating it as a hash set)
259
 * and recurses to its neighbors, if needed.
260
 *
261
 * @param  origin      Origin cell
262
 * @param  k           Maximum distance to move from the origin
263
 * @param  out         Array treated as a hash set, elements being either
264
 *                     H3Index or 0.
265
 * @param  distances   Scratch area, with elements paralleling the out array.
266
 *                     Elements indicate ijk distance from the origin cell to
267
 *                     the output cell
268
 * @param  maxIdx      Size of out and scratch arrays (must be
269
 * maxGridDiskSize(k))
270
 * @param  curK        Current distance from the origin
271
 */
272
H3Error _gridDiskDistancesInternal(H3Index origin, int k, H3Index *out,
273
156M
                                   int *distances, int64_t maxIdx, int curK) {
274
    // Put origin in the output array. out is used as a hash set.
275
156M
    int64_t off = origin % maxIdx;
276
4.63G
    while (out[off] != 0 && out[off] != origin) {
277
4.48G
        off = (off + 1) % maxIdx;
278
4.48G
    }
279
280
    // We either got a free slot in the hash set or hit a duplicate
281
    // We might need to process the duplicate anyways because we got
282
    // here on a longer path before.
283
156M
    if (out[off] == origin && distances[off] <= curK) return E_SUCCESS;
284
285
28.0M
    out[off] = origin;
286
28.0M
    distances[off] = curK;
287
288
    // Base case: reached an index k away from the origin.
289
28.0M
    if (curK >= k) return E_SUCCESS;
290
291
    // Recurse to all neighbors in no particular order.
292
182M
    for (int i = 0; i < 6; i++) {
293
156M
        int rotations = 0;
294
156M
        H3Index nextNeighbor;
295
156M
        H3Error neighborResult = h3NeighborRotations(origin, DIRECTIONS[i],
296
156M
                                                     &rotations, &nextNeighbor);
297
156M
        if (neighborResult != E_PENTAGON) {
298
            // E_PENTAGON is an expected case when trying to traverse off of
299
            // pentagons.
300
156M
            if (neighborResult != E_SUCCESS) {
301
455
                return neighborResult;
302
455
            }
303
156M
            neighborResult = _gridDiskDistancesInternal(
304
156M
                nextNeighbor, k, out, distances, maxIdx, curK + 1);
305
156M
            if (neighborResult) {
306
2.15k
                return neighborResult;
307
2.15k
            }
308
156M
        }
309
156M
    }
310
26.1M
    return E_SUCCESS;
311
26.1M
}
312
313
/**
314
 * Safe but slow version of gridDiskDistances (also called by it when needed).
315
 *
316
 * Adds the origin cell to the output set (treating it as a hash set)
317
 * and recurses to its neighbors, if needed.
318
 *
319
 * @param  origin      Origin cell
320
 * @param  k           Maximum distance to move from the origin
321
 * @param  out         Array treated as a hash set, elements being either
322
 *                     H3Index or 0.
323
 * @param  distances   Scratch area, with elements paralleling the out array.
324
 *                     Elements indicate ijk distance from the origin cell to
325
 *                     the output cell
326
 */
327
H3Error H3_EXPORT(gridDiskDistancesSafe)(H3Index origin, int k, H3Index *out,
328
371
                                         int *distances) {
329
371
    int64_t maxIdx;
330
371
    H3Error err = H3_EXPORT(maxGridDiskSize)(k, &maxIdx);
331
371
    if (err) {
332
0
        return err;
333
0
    }
334
371
    return _gridDiskDistancesInternal(origin, k, out, distances, maxIdx, 0);
335
371
}
336
337
/**
338
 * Maximum number of cells that result from the gridRing algorithm with the
339
 * given k.
340
 *
341
 * @param   k   k value, k >= 0.
342
 * @param out   size in indexes
343
 */
344
0
H3Error H3_EXPORT(maxGridRingSize)(int k, int64_t *out) {
345
0
    if (k < 0) {
346
0
        return E_DOMAIN;
347
0
    }
348
0
    if (k == 0) {
349
0
        *out = 1;
350
0
        return E_SUCCESS;
351
0
    }
352
0
    *out = 6 * (int64_t)k;
353
0
    return E_SUCCESS;
354
0
}
355
356
/**
357
 * Returns the "hollow" ring of hexagons at exactly grid distance k from
358
 * the origin hexagon. In particular, k=0 returns just the origin hexagon.
359
 *
360
 * Elements of the output array may be left zero, as can happen when crossing a
361
 * pentagon.
362
 *
363
 * @param origin Origin location.
364
 * @param k k >= 0
365
 * @param out Array which must be of size 6 * k (or 1 if k == 0)
366
 * @return 0 if successful; nonzero otherwise.
367
 */
368
371
H3Error H3_EXPORT(gridRing)(H3Index origin, int k, H3Index *out) {
369
    // Optimistically try the faster gridDiskUnsafe algorithm first
370
371
    const H3Error failed = H3_EXPORT(gridRingUnsafe)(origin, k, out);
371
371
    if (!failed) {
372
153
        return E_SUCCESS;
373
153
    }
374
    // Fast algo failed, fall back to slower, correct algo
375
    // and also wipe out array because contents untrustworthy
376
218
    memset(out, 0, 6 * k * sizeof(H3Index));
377
218
    return _gridRingInternal(origin, k, out);
378
371
}
379
380
/**
381
 * Internal algorithm for the safe but slow version of gridRing
382
 *
383
 * This function uses _gridDiskDistancesInternal to get all cells up to distance
384
 * k, then filters those results to include only cells exactly at distance k.
385
 *
386
 * @param origin Origin cell.
387
 * @param k Exact distance to filter for.
388
 * @param out Array which must be of size 6 * k (or 1 if k == 0). This is where
389
 * the final filtered results will be placed.
390
 * @return H3Error indicating success or failure.
391
 */
392
218
H3Error _gridRingInternal(H3Index origin, int k, H3Index *out) {
393
    // Short-circuit on 'identity' ring
394
218
    if (k == 0) {
395
0
        out[0] = origin;
396
0
        return E_SUCCESS;
397
0
    }
398
399
218
    int64_t maxIdx;
400
218
    H3Error err = H3_EXPORT(maxGridDiskSize)(k, &maxIdx);
401
218
    if (err) {
402
0
        return err;
403
0
    }
404
405
218
    H3Index *disk_out = H3_MEMORY(calloc)(maxIdx, sizeof(H3Index));
406
218
    if (!disk_out) {
407
0
        return E_MEMORY_ALLOC;
408
0
    }
409
218
    int *disk_distances = H3_MEMORY(calloc)(maxIdx, sizeof(int));
410
218
    if (!disk_distances) {
411
0
        H3_MEMORY(free)(disk_out);
412
0
        return E_MEMORY_ALLOC;
413
0
    }
414
415
218
    err = _gridDiskDistancesInternal(origin, k, disk_out, disk_distances,
416
218
                                     maxIdx, 0);
417
218
    if (err) {
418
87
        H3_MEMORY(free)(disk_out);
419
87
        H3_MEMORY(free)(disk_distances);
420
87
        return err;
421
87
    }
422
423
131
    int current_idx = 0;
424
468k
    for (int64_t i = 0; i < maxIdx; ++i) {
425
468k
        if (disk_out[i] != 0 && disk_distances[i] == k) {
426
12.9k
            out[current_idx++] = disk_out[i];
427
12.9k
        }
428
468k
    }
429
131
    H3_MEMORY(free)(disk_out);
430
131
    H3_MEMORY(free)(disk_distances);
431
131
    return E_SUCCESS;
432
218
}
433
434
/**
435
 * Returns the hexagon index neighboring the origin, in the direction dir.
436
 *
437
 * Implementation note: The only reachable case where this returns 0 is if the
438
 * origin is a pentagon and the translation is in the k direction. Thus,
439
 * 0 can only be returned if origin is a pentagon.
440
 *
441
 * @param origin Origin index
442
 * @param dir Direction to move in
443
 * @param rotations Number of ccw rotations to perform to reorient the
444
 *                  translation vector. Will be modified to the new number of
445
 *                  rotations to perform (such as when crossing a face edge.)
446
 * @param out H3Index of the specified neighbor if succesful
447
 * @return E_SUCCESS on success
448
 */
449
H3Error h3NeighborRotations(H3Index origin, Direction dir, int *rotations,
450
319M
                            H3Index *out) {
451
319M
    H3Index current = origin;
452
453
319M
    if (dir < CENTER_DIGIT || dir >= INVALID_DIGIT) {
454
602
        return E_FAILED;
455
602
    }
456
    // Ensure that rotations is modulo'd by 6 before any possible addition,
457
    // to protect against signed integer overflow.
458
319M
    *rotations = *rotations % 6;
459
329M
    for (int i = 0; i < *rotations; i++) {
460
9.43M
        dir = _rotate60ccw(dir);
461
9.43M
    }
462
463
319M
    int newRotations = 0;
464
319M
    int oldBaseCell = H3_GET_BASE_CELL(current);
465
319M
    if (NEVER(oldBaseCell < 0) || oldBaseCell >= NUM_BASE_CELLS) {
466
        // Base cells less than zero can not be represented in an index
467
542
        return E_CELL_INVALID;
468
542
    }
469
319M
    Direction oldLeadingDigit = _h3LeadingNonZeroDigit(current);
470
471
    // Adjust the indexing digits and, if needed, the base cell.
472
319M
    int r = H3_GET_RESOLUTION(current) - 1;
473
544M
    while (true) {
474
544M
        if (r == -1) {
475
20.3M
            H3_SET_BASE_CELL(current, baseCellNeighbors[oldBaseCell][dir]);
476
20.3M
            newRotations = baseCellNeighbor60CCWRots[oldBaseCell][dir];
477
478
20.3M
            if (H3_GET_BASE_CELL(current) == INVALID_BASE_CELL) {
479
                // Adjust for the deleted k vertex at the base cell level.
480
                // This edge actually borders a different neighbor.
481
117k
                H3_SET_BASE_CELL(current,
482
117k
                                 baseCellNeighbors[oldBaseCell][IK_AXES_DIGIT]);
483
117k
                newRotations =
484
117k
                    baseCellNeighbor60CCWRots[oldBaseCell][IK_AXES_DIGIT];
485
486
                // perform the adjustment for the k-subsequence we're skipping
487
                // over.
488
117k
                current = _h3Rotate60ccw(current);
489
117k
                *rotations = *rotations + 1;
490
117k
            }
491
492
20.3M
            break;
493
523M
        } else {
494
523M
            Direction oldDigit = H3_GET_INDEX_DIGIT(current, r + 1);
495
523M
            Direction nextDir;
496
523M
            if (oldDigit == INVALID_DIGIT) {
497
                // Only possible on invalid input
498
1.54k
                return E_CELL_INVALID;
499
523M
            } else if (isResolutionClassIII(r + 1)) {
500
256M
                H3_SET_INDEX_DIGIT(current, r + 1, NEW_DIGIT_II[oldDigit][dir]);
501
256M
                nextDir = NEW_ADJUSTMENT_II[oldDigit][dir];
502
267M
            } else {
503
267M
                H3_SET_INDEX_DIGIT(current, r + 1,
504
267M
                                   NEW_DIGIT_III[oldDigit][dir]);
505
267M
                nextDir = NEW_ADJUSTMENT_III[oldDigit][dir];
506
267M
            }
507
508
523M
            if (nextDir != CENTER_DIGIT) {
509
224M
                dir = nextDir;
510
224M
                r--;
511
299M
            } else {
512
                // No more adjustment to perform
513
299M
                break;
514
299M
            }
515
523M
        }
516
544M
    }
517
518
319M
    int newBaseCell = H3_GET_BASE_CELL(current);
519
319M
    if (_isBaseCellPentagon(newBaseCell)) {
520
59.4M
        int alreadyAdjustedKSubsequence = 0;
521
522
        // force rotation out of missing k-axes sub-sequence
523
59.4M
        if (_h3LeadingNonZeroDigit(current) == K_AXES_DIGIT) {
524
975k
            if (oldBaseCell != newBaseCell) {
525
                // in this case, we traversed into the deleted
526
                // k subsequence of a pentagon base cell.
527
                // We need to rotate out of that case depending
528
                // on how we got here.
529
                // check for a cw/ccw offset face; default is ccw
530
531
165k
                if (ALWAYS(_baseCellIsCwOffset(
532
165k
                        newBaseCell,
533
165k
                        baseCellData[oldBaseCell].homeFijk.face))) {
534
165k
                    current = _h3Rotate60cw(current);
535
165k
                } else {
536
                    // See cwOffsetPent in testGridDisk.c for why this is
537
                    // unreachable.
538
0
                    current = _h3Rotate60ccw(current);
539
0
                }
540
0
                alreadyAdjustedKSubsequence = 1;
541
810k
            } else {
542
                // In this case, we traversed into the deleted
543
                // k subsequence from within the same pentagon
544
                // base cell.
545
810k
                if (oldLeadingDigit == CENTER_DIGIT) {
546
                    // Undefined: the k direction is deleted from here
547
31.9k
                    return E_PENTAGON;
548
778k
                } else if (oldLeadingDigit == JK_AXES_DIGIT) {
549
                    // Rotate out of the deleted k subsequence
550
                    // We also need an additional change to the direction we're
551
                    // moving in
552
386k
                    current = _h3Rotate60ccw(current);
553
386k
                    *rotations = *rotations + 1;
554
392k
                } else if (oldLeadingDigit == IK_AXES_DIGIT) {
555
                    // Rotate out of the deleted k subsequence
556
                    // We also need an additional change to the direction we're
557
                    // moving in
558
391k
                    current = _h3Rotate60cw(current);
559
391k
                    *rotations = *rotations + 5;
560
391k
                } else {
561
                    // TODO: Should never occur, but is reachable by fuzzer
562
547
                    return E_FAILED;
563
547
                }
564
810k
            }
565
975k
        }
566
567
62.0M
        for (int i = 0; i < newRotations; i++)
568
2.55M
            current = _h3RotatePent60ccw(current);
569
570
        // Account for differing orientation of the base cells (this edge
571
        // might not follow properties of some other edges.)
572
59.4M
        if (oldBaseCell != newBaseCell) {
573
1.66M
            if (_isBaseCellPolarPentagon(newBaseCell)) {
574
                // 'polar' base cells behave differently because they have all
575
                // i neighbors.
576
403k
                if (oldBaseCell != 118 && oldBaseCell != 8 &&
577
403k
                    _h3LeadingNonZeroDigit(current) != JK_AXES_DIGIT) {
578
303k
                    *rotations = *rotations + 1;
579
303k
                }
580
1.26M
            } else if (_h3LeadingNonZeroDigit(current) == IK_AXES_DIGIT &&
581
1.26M
                       !alreadyAdjustedKSubsequence) {
582
                // account for distortion introduced to the 5 neighbor by the
583
                // deleted k subsequence.
584
153k
                *rotations = *rotations + 1;
585
153k
            }
586
1.66M
        }
587
260M
    } else {
588
289M
        for (int i = 0; i < newRotations; i++)
589
29.2M
            current = _h3Rotate60ccw(current);
590
260M
    }
591
592
319M
    *rotations = (*rotations + newRotations) % 6;
593
319M
    *out = current;
594
595
319M
    return E_SUCCESS;
596
319M
}
597
598
/**
599
 * Get the direction from the origin to a given neighbor. This is effectively
600
 * the reverse operation for h3NeighborRotations. Returns INVALID_DIGIT if the
601
 * cells are not neighbors.
602
 *
603
 * TODO: This is currently a brute-force algorithm, but as it's O(6) that's
604
 * probably acceptable.
605
 */
606
1.44k
Direction directionForNeighbor(H3Index origin, H3Index destination) {
607
1.44k
    bool isPent = H3_EXPORT(isPentagon)(origin);
608
    // Checks each neighbor, in order, to determine which direction the
609
    // destination neighbor is located. Skips CENTER_DIGIT since that
610
    // would be the origin; skips deleted K direction for pentagons.
611
1.44k
    for (Direction direction = isPent ? J_AXES_DIGIT : K_AXES_DIGIT;
612
9.50k
         direction < NUM_DIGITS; direction++) {
613
8.17k
        H3Index neighbor;
614
8.17k
        int rotations = 0;
615
8.17k
        H3Error neighborError =
616
8.17k
            h3NeighborRotations(origin, direction, &rotations, &neighbor);
617
8.17k
        if (!neighborError && neighbor == destination) {
618
110
            return direction;
619
110
        }
620
8.17k
    }
621
1.33k
    return INVALID_DIGIT;
622
1.44k
}
623
624
/**
625
 * gridDiskUnsafe produces indexes within k distance of the origin index.
626
 * Output behavior is undefined when one of the indexes returned by this
627
 * function is a pentagon or is in the pentagon distortion area.
628
 *
629
 * k-ring 0 is defined as the origin index, k-ring 1 is defined as k-ring 0 and
630
 * all neighboring indexes, and so on.
631
 *
632
 * Output is placed in the provided array in order of increasing distance from
633
 * the origin.
634
 *
635
 * @param origin Origin location.
636
 * @param k k >= 0
637
 * @param out Array which must be of size maxGridDiskSize(k).
638
 * @return 0 if no pentagon or pentagonal distortion area was encountered.
639
 */
640
62.3k
H3Error H3_EXPORT(gridDiskUnsafe)(H3Index origin, int k, H3Index *out) {
641
62.3k
    return H3_EXPORT(gridDiskDistancesUnsafe)(origin, k, out, NULL);
642
62.3k
}
643
644
/**
645
 * gridDiskDistancesUnsafe produces indexes within k distance of the origin
646
 * index. Output behavior is undefined when one of the indexes returned by this
647
 * function is a pentagon or is in the pentagon distortion area.
648
 *
649
 * k-ring 0 is defined as the origin index, k-ring 1 is defined as k-ring 0 and
650
 * all neighboring indexes, and so on.
651
 *
652
 * Output is placed in the provided array in order of increasing distance from
653
 * the origin. The distances in hexagons is placed in the distances array at
654
 * the same offset.
655
 *
656
 * @param origin Origin location.
657
 * @param k k >= 0
658
 * @param out Array which must be of size maxGridDiskSize(k).
659
 * @param distances Null or array which must be of size maxGridDiskSize(k).
660
 * @return 0 if no pentagon or pentagonal distortion area was encountered.
661
 */
662
H3Error H3_EXPORT(gridDiskDistancesUnsafe)(H3Index origin, int k, H3Index *out,
663
23.1M
                                           int *distances) {
664
    // Return codes:
665
    // 1 Pentagon was encountered
666
    // 2 Pentagon distortion (deleted k subsequence) was encountered
667
    // Pentagon being encountered is not itself a problem; really the deleted
668
    // k-subsequence is the problem, but for compatibility reasons we fail on
669
    // the pentagon.
670
23.1M
    if (k < 0) {
671
0
        return E_DOMAIN;
672
0
    }
673
674
    // k must be >= 0, so origin is always needed
675
23.1M
    int idx = 0;
676
23.1M
    out[idx] = origin;
677
23.1M
    if (distances) {
678
742
        distances[idx] = 0;
679
742
    }
680
23.1M
    idx++;
681
682
23.1M
    if (H3_EXPORT(isPentagon)(origin)) {
683
        // Pentagon was encountered; bail out as user doesn't want this.
684
6.38k
        return E_PENTAGON;
685
6.38k
    }
686
687
    // 0 < ring <= k, current ring
688
23.1M
    int ring = 1;
689
    // 0 <= direction < 6, current side of the ring
690
23.1M
    int direction = 0;
691
    // 0 <= i < ring, current position on the side of the ring
692
23.1M
    int i = 0;
693
    // Number of 60 degree ccw rotations to perform on the direction (based on
694
    // which faces have been crossed.)
695
23.1M
    int rotations = 0;
696
697
163M
    while (ring <= k) {
698
140M
        if (direction == 0 && i == 0) {
699
            // Not putting in the output set as it will be done later, at
700
            // the end of this ring.
701
23.1M
            H3Error neighborResult = h3NeighborRotations(
702
23.1M
                origin, NEXT_RING_DIRECTION, &rotations, &origin);
703
23.1M
            if (neighborResult) {
704
                // Should not be possible because `origin` would have to be a
705
                // pentagon
706
                // TODO: Reachable via fuzzer
707
307
                return neighborResult;
708
307
            }
709
710
23.1M
            if (H3_EXPORT(isPentagon)(origin)) {
711
                // Pentagon was encountered; bail out as user doesn't want this.
712
9.82k
                return E_PENTAGON;
713
9.82k
            }
714
23.1M
        }
715
716
140M
        H3Error neighborResult = h3NeighborRotations(
717
140M
            origin, DIRECTIONS[direction], &rotations, &origin);
718
140M
        if (neighborResult) {
719
84
            return neighborResult;
720
84
        }
721
140M
        out[idx] = origin;
722
140M
        if (distances) {
723
737k
            distances[idx] = ring;
724
737k
        }
725
140M
        idx++;
726
727
140M
        i++;
728
        // Check if end of this side of the k-ring
729
140M
        if (i == ring) {
730
138M
            i = 0;
731
138M
            direction++;
732
            // Check if end of this ring.
733
138M
            if (direction == 6) {
734
23.0M
                direction = 0;
735
23.0M
                ring++;
736
23.0M
            }
737
138M
        }
738
739
140M
        if (H3_EXPORT(isPentagon)(origin)) {
740
            // Pentagon was encountered; bail out as user doesn't want this.
741
24.2k
            return E_PENTAGON;
742
24.2k
        }
743
140M
    }
744
23.1M
    return E_SUCCESS;
745
23.1M
}
746
747
/**
748
 * gridDisksUnsafe takes an array of input hex IDs and a max k-ring and returns
749
 * an array of hexagon IDs sorted first by the original hex IDs and then by the
750
 * k-ring (0 to max), with no guaranteed sorting within each k-ring group.
751
 *
752
 * @param h3Set A pointer to an array of H3Indexes
753
 * @param length The total number of H3Indexes in h3Set
754
 * @param k The number of rings to generate
755
 * @param out A pointer to the output memory to dump the new set of H3Indexes to
756
 *            The memory block should be equal to maxGridDiskSize(k) * length
757
 * @return 0 if no pentagon is encountered. Cannot trust output otherwise
758
 */
759
H3Error H3_EXPORT(gridDisksUnsafe)(H3Index *h3Set, int length, int k,
760
349
                                   H3Index *out) {
761
349
    H3Index *segment;
762
349
    int64_t segmentSize;
763
349
    H3Error err = H3_EXPORT(maxGridDiskSize)(k, &segmentSize);
764
349
    if (err) {
765
0
        return err;
766
0
    }
767
62.3k
    for (int i = 0; i < length; i++) {
768
        // Determine the appropriate segment of the output array to operate on
769
62.0k
        segment = out + i * segmentSize;
770
62.0k
        H3Error failed = H3_EXPORT(gridDiskUnsafe)(h3Set[i], k, segment);
771
62.0k
        if (failed) return failed;
772
62.0k
    }
773
296
    return E_SUCCESS;
774
349
}
775
776
/**
777
 * Returns the "hollow" ring of hexagons at exactly grid distance k from
778
 * the origin hexagon. In particular, k=0 returns just the origin hexagon.
779
 *
780
 * A nonzero failure code may be returned in some cases, for example,
781
 * if a pentagon is encountered.
782
 * Failure cases may be fixed in future versions.
783
 *
784
 * @param origin Origin location.
785
 * @param k k >= 0
786
 * @param out Array which must be of size 6 * k (or 1 if k == 0)
787
 * @return 0 if successful; nonzero otherwise.
788
 */
789
742
H3Error H3_EXPORT(gridRingUnsafe)(H3Index origin, int k, H3Index *out) {
790
742
    if (k < 0) {
791
0
        return E_DOMAIN;
792
0
    }
793
    // Short-circuit on 'identity' ring
794
742
    if (k == 0) {
795
122
        out[0] = origin;
796
122
        return E_SUCCESS;
797
122
    }
798
620
    int idx = 0;
799
    // Number of 60 degree ccw rotations to perform on the direction (based on
800
    // which faces have been crossed.)
801
620
    int rotations = 0;
802
    // Scratch structure for checking for pentagons
803
620
    if (H3_EXPORT(isPentagon)(origin)) {
804
        // Pentagon was encountered; bail out as user doesn't want this.
805
18
        return E_PENTAGON;
806
18
    }
807
808
13.1k
    for (int ring = 0; ring < k; ring++) {
809
12.7k
        H3Error neighborResult = h3NeighborRotations(
810
12.7k
            origin, NEXT_RING_DIRECTION, &rotations, &origin);
811
12.7k
        if (neighborResult) {
812
            // Should not be possible because `origin` would have to be a
813
            // pentagon
814
            // TODO: Reachable via fuzzer
815
130
            return neighborResult;
816
130
        }
817
818
12.6k
        if (H3_EXPORT(isPentagon)(origin)) {
819
54
            return E_PENTAGON;
820
54
        }
821
12.6k
    }
822
823
418
    H3Index lastIndex = origin;
824
825
418
    out[idx] = origin;
826
418
    idx++;
827
828
2.58k
    for (int direction = 0; direction < 6; direction++) {
829
66.5k
        for (int pos = 0; pos < k; pos++) {
830
64.4k
            H3Error neighborResult = h3NeighborRotations(
831
64.4k
                origin, DIRECTIONS[direction], &rotations, &origin);
832
64.4k
            if (neighborResult) {
833
                // Should not be possible because `origin` would have to be a
834
                // pentagon
835
                // TODO: Reachable via fuzzer
836
18
                return neighborResult;
837
18
            }
838
839
            // Skip the very last index, it was already added. We do
840
            // however need to traverse to it because of the pentagonal
841
            // distortion check, below.
842
64.3k
            if (pos != k - 1 || direction != 5) {
843
64.0k
                out[idx] = origin;
844
64.0k
                idx++;
845
846
64.0k
                if (H3_EXPORT(isPentagon)(origin)) {
847
58
                    return E_PENTAGON;
848
58
                }
849
64.0k
            }
850
64.3k
        }
851
2.24k
    }
852
853
    // Check that this matches the expected lastIndex, if it doesn't,
854
    // it indicates pentagonal distortion occurred and we should report
855
    // failure.
856
342
    if (lastIndex != origin) {
857
158
        return E_PENTAGON;
858
158
    }
859
184
    return E_SUCCESS;
860
342
}
861
862
/**
863
 * maxPolygonToCellsSize returns the number of cells to allocate space for
864
 * when performing a polygonToCells on the given GeoJSON-like data structure.
865
 *
866
 * The size is the maximum of either the number of points in the geoloop or the
867
 * number of cells in the bounding box of the geoloop.
868
 *
869
 * @param geoPolygon A GeoJSON-like data structure indicating the poly to fill
870
 * @param res Hexagon resolution (0-15)
871
 * @param out number of cells to allocate for
872
 * @return 0 (E_SUCCESS) on success.
873
 */
874
H3Error H3_EXPORT(maxPolygonToCellsSize)(const GeoPolygon *geoPolygon, int res,
875
9.58k
                                         uint32_t flags, int64_t *out) {
876
9.58k
    H3Error flagErr = validatePolygonFlags(flags);
877
9.58k
    if (flagErr) {
878
0
        return flagErr;
879
0
    }
880
    // Get the bounding box for the GeoJSON-like struct
881
9.58k
    BBox bbox;
882
9.58k
    const GeoLoop geoloop = geoPolygon->geoloop;
883
9.58k
    bboxFromGeoLoop(&geoloop, &bbox);
884
9.58k
    int64_t numHexagons;
885
9.58k
    H3Error estimateErr = bboxHexEstimate(&bbox, res, &numHexagons);
886
9.58k
    if (estimateErr) {
887
232
        return estimateErr;
888
232
    }
889
    // This algorithm assumes that the number of vertices is usually less than
890
    // the number of hexagons, but when it's wrong, this will keep it from
891
    // failing
892
9.35k
    int totalVerts = geoloop.numVerts;
893
91.9k
    for (int i = 0; i < geoPolygon->numHoles; i++) {
894
82.5k
        totalVerts += geoPolygon->holes[i].numVerts;
895
82.5k
    }
896
9.35k
    if (numHexagons < totalVerts) numHexagons = totalVerts;
897
    // When the polygon is very small, near an icosahedron edge and is an odd
898
    // resolution, the line tracing needs an extra buffer than the estimator
899
    // function provides (but beefing that up to cover causes most situations to
900
    // overallocate memory)
901
9.35k
    numHexagons += POLYGON_TO_CELLS_BUFFER;
902
9.35k
    *out = numHexagons;
903
9.35k
    return E_SUCCESS;
904
9.58k
}
905
906
/**
907
 * _getEdgeHexagons takes a given geoloop ring (either the main geoloop or
908
 * one of the holes) and traces it with hexagons and updates the search and
909
 * found memory blocks. This is used for determining the initial hexagon set
910
 * for the polygonToCells algorithm to execute on.
911
 *
912
 * @param geoloop The geoloop (or hole) to be traced
913
 * @param numHexagons The maximum number of hexagons possible for the geoloop
914
 *                    (also the bounds of the search and found arrays)
915
 * @param res The hexagon resolution (0-15)
916
 * @param numSearchHexes The number of hexagons found so far to be searched
917
 * @param search The block of memory containing the hexagons to search from
918
 * @param found The block of memory containing the hexagons found from the
919
 * search
920
 *
921
 * @return An error code if the hash function cannot insert a found hexagon
922
 *         into the found array.
923
 */
924
H3Error _getEdgeHexagons(const GeoLoop *geoloop, int64_t numHexagons, int res,
925
                         int64_t *numSearchHexes, H3Index *search,
926
31.9k
                         H3Index *found) {
927
2.62M
    for (int i = 0; i < geoloop->numVerts; i++) {
928
2.59M
        LatLng origin = geoloop->verts[i];
929
2.59M
        LatLng destination = i == geoloop->numVerts - 1 ? geoloop->verts[0]
930
2.59M
                                                        : geoloop->verts[i + 1];
931
2.59M
        int64_t numHexesEstimate;
932
2.59M
        H3Error estimateErr =
933
2.59M
            lineHexEstimate(&origin, &destination, res, &numHexesEstimate);
934
2.59M
        if (estimateErr) {
935
472
            return estimateErr;
936
472
        }
937
56.5M
        for (int64_t j = 0; j < numHexesEstimate; j++) {
938
53.9M
            LatLng interpolate;
939
53.9M
            double invNumHexesEst = 1.0 / numHexesEstimate;
940
53.9M
            interpolate.lat =
941
53.9M
                (origin.lat * (numHexesEstimate - j) * invNumHexesEst) +
942
53.9M
                (destination.lat * j * invNumHexesEst);
943
53.9M
            interpolate.lng =
944
53.9M
                (origin.lng * (numHexesEstimate - j) * invNumHexesEst) +
945
53.9M
                (destination.lng * j * invNumHexesEst);
946
53.9M
            H3Index pointHex;
947
53.9M
            H3Error e = H3_EXPORT(latLngToCell)(&interpolate, res, &pointHex);
948
53.9M
            if (e) {
949
144
                return e;
950
144
            }
951
            // A simple hash to store the hexagon, or move to another place if
952
            // needed
953
53.9M
            int64_t loc = (int64_t)(pointHex % numHexagons);
954
53.9M
            int64_t loopCount = 0;
955
1.98G
            while (found[loc] != 0) {
956
                // If this conditional is reached, the `found` memory block is
957
                // too small for the given polygon. This should not happen.
958
                // TODO: Reachable via fuzzer
959
1.95G
                if (loopCount > numHexagons) return E_FAILED;
960
1.95G
                if (found[loc] == pointHex)
961
27.1M
                    break;  // At least two points of the geoloop index to the
962
                            // same cell
963
1.93G
                loc = (loc + 1) % numHexagons;
964
1.93G
                loopCount++;
965
1.93G
            }
966
53.9M
            if (found[loc] == pointHex)
967
27.1M
                continue;  // Skip this hex, already exists in the found hash
968
            // Otherwise, set it in the found hash for now
969
26.7M
            found[loc] = pointHex;
970
971
26.7M
            search[*numSearchHexes] = pointHex;
972
26.7M
            (*numSearchHexes)++;
973
26.7M
        }
974
2.59M
    }
975
30.9k
    return E_SUCCESS;
976
31.9k
}
977
978
/**
979
 * polygonToCells takes a given GeoJSON-like data structure and preallocated,
980
 * zeroed memory, and fills it with the hexagons that are contained by
981
 * the GeoJSON-like data structure.
982
 *
983
 * This implementation traces the GeoJSON geoloop(s) in cartesian space with
984
 * hexagons, tests them and their neighbors to be contained by the geoloop(s),
985
 * and then any newly found hexagons are used to test again until no new
986
 * hexagons are found.
987
 *
988
 * @param geoPolygon The geoloop and holes defining the relevant area
989
 * @param res The Hexagon resolution (0-15)
990
 * @param out The slab of zeroed memory to write to. Assumed to be big enough.
991
 */
992
H3Error H3_EXPORT(polygonToCells)(const GeoPolygon *geoPolygon, int res,
993
4.33k
                                  uint32_t flags, H3Index *out) {
994
4.33k
    H3Error flagErr = validatePolygonFlags(flags);
995
4.33k
    if (flagErr) {
996
0
        return flagErr;
997
0
    }
998
    // One of the goals of the polygonToCells algorithm is that two adjacent
999
    // polygons with zero overlap have zero overlapping hexagons. That the
1000
    // hexagons are uniquely assigned. There are a few approaches to take here,
1001
    // such as deciding based on which polygon has the greatest overlapping area
1002
    // of the hexagon, or the most number of contained points on the hexagon
1003
    // (using the center point as a tiebreaker).
1004
    //
1005
    // But if the polygons are convex, both of these more complex algorithms can
1006
    // be reduced down to checking whether or not the center of the hexagon is
1007
    // contained in the polygon, and so this is the approach that this
1008
    // polygonToCells algorithm will follow, as it's simpler, faster, and the
1009
    // error for concave polygons is still minimal (only affecting concave
1010
    // shapes on the order of magnitude of the hexagon size or smaller, not
1011
    // impacting larger concave shapes)
1012
    //
1013
    // This first part is identical to the maxPolygonToCellsSize above.
1014
1015
    // Get the bounding boxes for the polygon and any holes
1016
4.33k
    BBox *bboxes = H3_MEMORY(malloc)((geoPolygon->numHoles + 1) * sizeof(BBox));
1017
4.33k
    if (!bboxes) {
1018
0
        return E_MEMORY_ALLOC;
1019
0
    }
1020
4.33k
    bboxesFromGeoPolygon(geoPolygon, bboxes);
1021
1022
    // Get the estimated number of hexagons and allocate some temporary memory
1023
    // for the hexagons
1024
4.33k
    int64_t numHexagons;
1025
4.33k
    H3Error numHexagonsError =
1026
4.33k
        H3_EXPORT(maxPolygonToCellsSize)(geoPolygon, res, flags, &numHexagons);
1027
4.33k
    if (numHexagonsError) {
1028
0
        H3_MEMORY(free)(bboxes);
1029
0
        return numHexagonsError;
1030
0
    }
1031
4.33k
    H3Index *search = H3_MEMORY(calloc)(numHexagons, sizeof(H3Index));
1032
4.33k
    if (!search) {
1033
0
        H3_MEMORY(free)(bboxes);
1034
0
        return E_MEMORY_ALLOC;
1035
0
    }
1036
4.33k
    H3Index *found = H3_MEMORY(calloc)(numHexagons, sizeof(H3Index));
1037
4.33k
    if (!found) {
1038
0
        H3_MEMORY(free)(bboxes);
1039
0
        H3_MEMORY(free)(search);
1040
0
        return E_MEMORY_ALLOC;
1041
0
    }
1042
1043
    // Some metadata for tracking the state of the search and found memory
1044
    // blocks
1045
4.33k
    int64_t numSearchHexes = 0;
1046
4.33k
    int64_t numFoundHexes = 0;
1047
1048
    // 1. Trace the hexagons along the polygon defining the outer geoloop and
1049
    // add them to the search hash. The hexagon containing the geoloop point
1050
    // may or may not be contained by the geoloop (as the hexagon's center
1051
    // point may be outside of the boundary.)
1052
4.33k
    const GeoLoop geoloop = geoPolygon->geoloop;
1053
4.33k
    H3Error edgeHexError = _getEdgeHexagons(&geoloop, numHexagons, res,
1054
4.33k
                                            &numSearchHexes, search, found);
1055
    // If this branch is reached, we have exceeded the maximum number of
1056
    // hexagons possible and need to clean up the allocated memory.
1057
    // TODO: Reachable via fuzzer
1058
4.33k
    if (edgeHexError) {
1059
900
        H3_MEMORY(free)(search);
1060
900
        H3_MEMORY(free)(found);
1061
900
        H3_MEMORY(free)(bboxes);
1062
900
        return edgeHexError;
1063
900
    }
1064
1065
    // 2. Iterate over all holes, trace the polygons defining the holes with
1066
    // hexagons and add to only the search hash. We're going to temporarily use
1067
    // the `found` hash to use for dedupe purposes and then re-zero it once
1068
    // we're done here, otherwise we'd have to scan the whole set on each insert
1069
    // to make sure there's no duplicates, which is very inefficient.
1070
30.9k
    for (int i = 0; i < geoPolygon->numHoles; i++) {
1071
27.6k
        GeoLoop *hole = &(geoPolygon->holes[i]);
1072
27.6k
        edgeHexError = _getEdgeHexagons(hole, numHexagons, res, &numSearchHexes,
1073
27.6k
                                        search, found);
1074
        // If this branch is reached, we have exceeded the maximum number of
1075
        // hexagons possible and need to clean up the allocated memory.
1076
        // TODO: Reachable via fuzzer
1077
27.6k
        if (edgeHexError) {
1078
156
            H3_MEMORY(free)(search);
1079
156
            H3_MEMORY(free)(found);
1080
156
            H3_MEMORY(free)(bboxes);
1081
156
            return edgeHexError;
1082
156
        }
1083
27.6k
    }
1084
1085
    // 3. Re-zero the found hash so it can be used in the main loop below
1086
1.73G
    for (int64_t i = 0; i < numHexagons; i++) found[i] = H3_NULL;
1087
1088
    // 4. Begin main loop. While the search hash is not empty do the following
1089
23.4k
    while (numSearchHexes > 0) {
1090
        // Iterate through all hexagons in the current search hash, then loop
1091
        // through all neighbors and test Point-in-Poly, if point-in-poly
1092
        // succeeds, add to out and found hashes if not already there.
1093
20.8k
        int64_t currentSearchNum = 0;
1094
20.8k
        int64_t i = 0;
1095
23.1M
        while (currentSearchNum < numSearchHexes) {
1096
23.1M
            H3Index ring[MAX_ONE_RING_SIZE] = {0};
1097
23.1M
            H3Index searchHex = search[i];
1098
23.1M
            H3_EXPORT(gridDisk)(searchHex, 1, ring);
1099
184M
            for (int j = 0; j < MAX_ONE_RING_SIZE; j++) {
1100
161M
                if (ring[j] == H3_NULL) {
1101
6.32k
                    continue;  // Skip if this was a pentagon and only had 5
1102
                               // neighbors
1103
6.32k
                }
1104
1105
161M
                H3Index hex = ring[j];
1106
1107
                // A simple hash to store the hexagon, or move to another place
1108
                // if needed. This MUST be done before the point-in-poly check
1109
                // since that's far more expensive
1110
161M
                int64_t loc = (int64_t)(hex % numHexagons);
1111
161M
                int64_t loopCount = 0;
1112
1.59G
                while (out[loc] != 0) {
1113
                    // If this branch is reached, we have exceeded the maximum
1114
                    // number of hexagons possible and need to clean up the
1115
                    // allocated memory.
1116
                    // TODO: Reachable via fuzzer
1117
1.51G
                    if (loopCount > numHexagons) {
1118
604
                        H3_MEMORY(free)(search);
1119
604
                        H3_MEMORY(free)(found);
1120
604
                        H3_MEMORY(free)(bboxes);
1121
604
                        return E_FAILED;
1122
604
                    }
1123
1.51G
                    if (out[loc] == hex) break;  // Skip duplicates found
1124
1.43G
                    loc = (loc + 1) % numHexagons;
1125
1.43G
                    loopCount++;
1126
1.43G
                }
1127
161M
                if (out[loc] == hex) {
1128
80.8M
                    continue;  // Skip this hex, already exists in the out hash
1129
80.8M
                }
1130
1131
                // Check if the hexagon is in the polygon or not
1132
80.9M
                LatLng hexCenter;
1133
80.9M
                H3_EXPORT(cellToLatLng)(hex, &hexCenter);
1134
1135
                // If not, skip
1136
80.9M
                if (!pointInsidePolygon(geoPolygon, bboxes, &hexCenter)) {
1137
66.4M
                    continue;
1138
66.4M
                }
1139
1140
                // Otherwise set it in the output array
1141
14.4M
                out[loc] = hex;
1142
1143
                // Set the hexagon in the found hash
1144
14.4M
                found[numFoundHexes] = hex;
1145
14.4M
                numFoundHexes++;
1146
14.4M
            }
1147
23.1M
            currentSearchNum++;
1148
23.1M
            i++;
1149
23.1M
        }
1150
1151
        // Swap the search and found pointers, copy the found hex count to the
1152
        // search hex count, and zero everything related to the found memory.
1153
20.2k
        H3Index *temp = search;
1154
20.2k
        search = found;
1155
20.2k
        found = temp;
1156
22.4M
        for (int64_t j = 0; j < numSearchHexes; j++) found[j] = 0;
1157
20.2k
        numSearchHexes = numFoundHexes;
1158
20.2k
        numFoundHexes = 0;
1159
        // Repeat until no new hexagons are found
1160
20.2k
    }
1161
    // The out memory structure should be complete, end it here
1162
2.67k
    H3_MEMORY(free)(bboxes);
1163
2.67k
    H3_MEMORY(free)(search);
1164
2.67k
    H3_MEMORY(free)(found);
1165
2.67k
    return E_SUCCESS;
1166
3.28k
}
1167
1168
/**
1169
 * Internal: Create a vertex graph from a set of hexagons. It is the
1170
 * responsibility of the caller to call destroyVertexGraph on the populated
1171
 * graph, otherwise the memory in the graph nodes will not be freed.
1172
 * @private
1173
 * @param h3Set    Set of hexagons
1174
 * @param numHexes Number of hexagons in the set
1175
 * @param graph    Output graph
1176
 */
1177
H3Error h3SetToVertexGraph(const H3Index *h3Set, const int numHexes,
1178
1.58k
                           VertexGraph *graph) {
1179
1.58k
    CellBoundary vertices;
1180
1.58k
    LatLng *fromVtx;
1181
1.58k
    LatLng *toVtx;
1182
1.58k
    VertexNode *edge;
1183
1.58k
    if (numHexes < 1) {
1184
        // We still need to init the graph, or calls to destroyVertexGraph will
1185
        // fail
1186
1
        initVertexGraph(graph, 0, 0);
1187
1
        return E_SUCCESS;
1188
1
    }
1189
1.58k
    int res = H3_GET_RESOLUTION(h3Set[0]);
1190
1.58k
    const int minBuckets = 6;
1191
    // TODO: Better way to calculate/guess?
1192
1.58k
    int numBuckets = numHexes > minBuckets ? numHexes : minBuckets;
1193
1.58k
    initVertexGraph(graph, numBuckets, res);
1194
    // Iterate through every hexagon
1195
1.90M
    for (int i = 0; i < numHexes; i++) {
1196
1.90M
        H3Error boundaryErr = H3_EXPORT(cellToBoundary)(h3Set[i], &vertices);
1197
1.90M
        if (boundaryErr) {
1198
            // Destroy vertex graph as caller will not know to do so.
1199
164
            destroyVertexGraph(graph);
1200
164
            return boundaryErr;
1201
164
        }
1202
        // iterate through every edge
1203
13.3M
        for (int j = 0; j < vertices.numVerts; j++) {
1204
11.4M
            fromVtx = &vertices.verts[j];
1205
11.4M
            toVtx = &vertices.verts[(j + 1) % vertices.numVerts];
1206
            // If we've seen this edge already, it will be reversed
1207
11.4M
            edge = findNodeForEdge(graph, toVtx, fromVtx);
1208
11.4M
            if (edge != NULL) {
1209
                // If we've seen it, drop it. No edge is shared by more than 2
1210
                // hexagons, so we'll never see it again.
1211
65.8k
                removeVertexNode(graph, edge);
1212
11.4M
            } else {
1213
                // Add a new node for this edge
1214
11.4M
                addVertexNode(graph, fromVtx, toVtx);
1215
11.4M
            }
1216
11.4M
        }
1217
1.90M
    }
1218
1.41k
    return E_SUCCESS;
1219
1.58k
}
1220
1221
/**
1222
 * Internal: Create a LinkedGeoPolygon from a vertex graph. It is the
1223
 * responsibility of the caller to call destroyLinkedMultiPolygon on the
1224
 * populated linked geo structure, or the memory for that structure will not be
1225
 * freed.
1226
 * @private
1227
 * @param graph Input graph
1228
 * @param out   Output polygon
1229
 */
1230
1.42k
void _vertexGraphToLinkedGeo(VertexGraph *graph, LinkedGeoPolygon *out) {
1231
1.42k
    *out = (LinkedGeoPolygon){0};
1232
1.42k
    LinkedGeoLoop *loop;
1233
1.42k
    VertexNode *edge;
1234
1.42k
    LatLng nextVtx;
1235
    // Find the next unused entry point
1236
64.1k
    while ((edge = firstVertexNode(graph)) != NULL) {
1237
62.7k
        loop = addNewLinkedLoop(out);
1238
        // Walk the graph to get the outline
1239
298k
        do {
1240
298k
            addLinkedCoord(loop, &edge->from);
1241
298k
            nextVtx = edge->to;
1242
            // Remove frees the node, so we can't use edge after this
1243
298k
            removeVertexNode(graph, edge);
1244
298k
            edge = findNodeForVertex(graph, &nextVtx);
1245
298k
        } while (edge);
1246
62.7k
    }
1247
1.42k
}
1248
1249
/**
1250
 * Create a LinkedGeoPolygon describing the outline(s) of a set of  hexagons.
1251
 * Polygon outlines will follow GeoJSON MultiPolygon order: Each polygon will
1252
 * have one outer loop, which is first in the list, followed by any holes.
1253
 *
1254
 * It is the responsibility of the caller to call destroyLinkedMultiPolygon on
1255
 * the populated linked geo structure, or the memory for that structure will not
1256
 * be freed.
1257
 *
1258
 * It is expected that all hexagons in the set have the same resolution and
1259
 * that the set contains no duplicates. Behavior is undefined if duplicates
1260
 * or multiple resolutions are present, and the algorithm may produce
1261
 * unexpected or invalid output.
1262
 *
1263
 * @param h3Set    Set of hexagons
1264
 * @param numHexes Number of hexagons in set
1265
 * @param out      Output polygon
1266
 */
1267
H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set,
1268
                                             const int numHexes,
1269
882
                                             LinkedGeoPolygon *out) {
1270
882
    VertexGraph graph;
1271
882
    H3Error err = h3SetToVertexGraph(h3Set, numHexes, &graph);
1272
882
    if (err) {
1273
48
        return err;
1274
48
    }
1275
834
    _vertexGraphToLinkedGeo(&graph, out);
1276
834
    destroyVertexGraph(&graph);
1277
834
    H3Error normalizeResult = normalizeMultiPolygon(out);
1278
834
    if (normalizeResult) {
1279
170
        H3_EXPORT(destroyLinkedMultiPolygon)(out);
1280
170
    }
1281
834
    return normalizeResult;
1282
882
}