1/*
2 * Copyright 2016-2018 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#include <float.h>
22#include <math.h>
23#include <stdbool.h>
24#include <stdlib.h>
25#include "baseCells.h"
26#include "bbox.h"
27#include "faceijk.h"
28#include "geoCoord.h"
29#include "h3Index.h"
30#include "h3api.h"
31#include "linkedGeo.h"
32#include "polygon.h"
33#include "stackAlloc.h"
34#include "vertexGraph.h"
35
36/*
37 * Return codes from hexRange and related functions.
38 */
39
40#define HEX_RANGE_SUCCESS 0
41#define HEX_RANGE_PENTAGON 1
42#define HEX_RANGE_K_SUBSEQUENCE 2
43
44/**
45 * Directions used for traversing a hexagonal ring counterclockwise around
46 * {1, 0, 0}
47 *
48 * <pre>
49 * _
50 * _/ \\_
51 * / \\5/ \\
52 * \\0/ \\4/
53 * / \\_/ \\
54 * \\1/ \\3/
55 * \\2/
56 * </pre>
57 */
58static const Direction DIRECTIONS[6] = {J_AXES_DIGIT, JK_AXES_DIGIT,
59 K_AXES_DIGIT, IK_AXES_DIGIT,
60 I_AXES_DIGIT, IJ_AXES_DIGIT};
61
62/**
63 * Direction used for traversing to the next outward hexagonal ring.
64 */
65static const Direction NEXT_RING_DIRECTION = I_AXES_DIGIT;
66
67/**
68 * New digit when traversing along class II grids.
69 *
70 * Current digit -> direction -> new digit.
71 */
72static const Direction NEW_DIGIT_II[7][7] = {
73 {CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT,
74 IK_AXES_DIGIT, IJ_AXES_DIGIT},
75 {K_AXES_DIGIT, I_AXES_DIGIT, JK_AXES_DIGIT, IJ_AXES_DIGIT, IK_AXES_DIGIT,
76 J_AXES_DIGIT, CENTER_DIGIT},
77 {J_AXES_DIGIT, JK_AXES_DIGIT, K_AXES_DIGIT, I_AXES_DIGIT, IJ_AXES_DIGIT,
78 CENTER_DIGIT, IK_AXES_DIGIT},
79 {JK_AXES_DIGIT, IJ_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT, CENTER_DIGIT,
80 K_AXES_DIGIT, J_AXES_DIGIT},
81 {I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, J_AXES_DIGIT,
82 JK_AXES_DIGIT, K_AXES_DIGIT},
83 {IK_AXES_DIGIT, J_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, JK_AXES_DIGIT,
84 IJ_AXES_DIGIT, I_AXES_DIGIT},
85 {IJ_AXES_DIGIT, CENTER_DIGIT, IK_AXES_DIGIT, J_AXES_DIGIT, K_AXES_DIGIT,
86 I_AXES_DIGIT, JK_AXES_DIGIT}};
87
88/**
89 * New traversal direction when traversing along class II grids.
90 *
91 * Current digit -> direction -> new ap7 move (at coarser level).
92 */
93static const Direction NEW_ADJUSTMENT_II[7][7] = {
94 {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT,
95 CENTER_DIGIT, CENTER_DIGIT},
96 {CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT,
97 IK_AXES_DIGIT, CENTER_DIGIT},
98 {CENTER_DIGIT, CENTER_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
99 CENTER_DIGIT, J_AXES_DIGIT},
100 {CENTER_DIGIT, K_AXES_DIGIT, JK_AXES_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
101 CENTER_DIGIT, CENTER_DIGIT},
102 {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
103 I_AXES_DIGIT, IJ_AXES_DIGIT},
104 {CENTER_DIGIT, IK_AXES_DIGIT, CENTER_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
105 IK_AXES_DIGIT, CENTER_DIGIT},
106 {CENTER_DIGIT, CENTER_DIGIT, J_AXES_DIGIT, CENTER_DIGIT, IJ_AXES_DIGIT,
107 CENTER_DIGIT, IJ_AXES_DIGIT}};
108
109/**
110 * New traversal direction when traversing along class III grids.
111 *
112 * Current digit -> direction -> new ap7 move (at coarser level).
113 */
114static const Direction NEW_DIGIT_III[7][7] = {
115 {CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT,
116 IK_AXES_DIGIT, IJ_AXES_DIGIT},
117 {K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT,
118 IJ_AXES_DIGIT, CENTER_DIGIT},
119 {J_AXES_DIGIT, JK_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT,
120 CENTER_DIGIT, K_AXES_DIGIT},
121 {JK_AXES_DIGIT, I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT,
122 K_AXES_DIGIT, J_AXES_DIGIT},
123 {I_AXES_DIGIT, IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT,
124 J_AXES_DIGIT, JK_AXES_DIGIT},
125 {IK_AXES_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT,
126 JK_AXES_DIGIT, I_AXES_DIGIT},
127 {IJ_AXES_DIGIT, CENTER_DIGIT, K_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT,
128 I_AXES_DIGIT, IK_AXES_DIGIT}};
129
130/**
131 * New traversal direction when traversing along class III grids.
132 *
133 * Current digit -> direction -> new ap7 move (at coarser level).
134 */
135static const Direction NEW_ADJUSTMENT_III[7][7] = {
136 {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT,
137 CENTER_DIGIT, CENTER_DIGIT},
138 {CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
139 K_AXES_DIGIT, CENTER_DIGIT},
140 {CENTER_DIGIT, CENTER_DIGIT, J_AXES_DIGIT, J_AXES_DIGIT, CENTER_DIGIT,
141 CENTER_DIGIT, IJ_AXES_DIGIT},
142 {CENTER_DIGIT, JK_AXES_DIGIT, J_AXES_DIGIT, JK_AXES_DIGIT, CENTER_DIGIT,
143 CENTER_DIGIT, CENTER_DIGIT},
144 {CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
145 IK_AXES_DIGIT, I_AXES_DIGIT},
146 {CENTER_DIGIT, K_AXES_DIGIT, CENTER_DIGIT, CENTER_DIGIT, IK_AXES_DIGIT,
147 IK_AXES_DIGIT, CENTER_DIGIT},
148 {CENTER_DIGIT, CENTER_DIGIT, IJ_AXES_DIGIT, CENTER_DIGIT, I_AXES_DIGIT,
149 CENTER_DIGIT, IJ_AXES_DIGIT}};
150
151/**
152 * Maximum number of indices that result from the kRing algorithm with the given
153 * k. Formula source and proof: https://oeis.org/A003215
154 *
155 * @param k k value, k >= 0.
156 */
157int H3_EXPORT(maxKringSize)(int k) { return 3 * k * (k + 1) + 1; }
158
159/**
160 * k-rings produces indices within k distance of the origin index.
161 *
162 * k-ring 0 is defined as the origin index, k-ring 1 is defined as k-ring 0 and
163 * all neighboring indices, and so on.
164 *
165 * Output is placed in the provided array in no particular order. Elements of
166 * the output array may be left zero, as can happen when crossing a pentagon.
167 *
168 * @param origin Origin location.
169 * @param k k >= 0
170 * @param out Zero-filled array which must be of size maxKringSize(k).
171 */
172void H3_EXPORT(kRing)(H3Index origin, int k, H3Index* out) {
173 int maxIdx = H3_EXPORT(maxKringSize)(k);
174 int* distances = malloc(maxIdx * sizeof(int));
175 H3_EXPORT(kRingDistances)(origin, k, out, distances);
176 free(distances);
177}
178
179/**
180 * k-rings produces indices within k distance of the origin index.
181 *
182 * k-ring 0 is defined as the origin index, k-ring 1 is defined as k-ring 0 and
183 * all neighboring indices, and so on.
184 *
185 * Output is placed in the provided array in no particular order. Elements of
186 * the output array may be left zero, as can happen when crossing a pentagon.
187 *
188 * @param origin Origin location.
189 * @param k k >= 0
190 * @param out Zero-filled array which must be of size maxKringSize(k).
191 * @param distances Zero-filled array which must be of size maxKringSize(k).
192 */
193void H3_EXPORT(kRingDistances)(H3Index origin, int k, H3Index* out,
194 int* distances) {
195 const int maxIdx = H3_EXPORT(maxKringSize)(k);
196 // Optimistically try the faster hexRange algorithm first
197 const bool failed = H3_EXPORT(hexRangeDistances)(origin, k, out, distances);
198 if (failed) {
199 // Fast algo failed, fall back to slower, correct algo
200 // and also wipe out array because contents untrustworthy
201 memset(out, 0, maxIdx * sizeof(out[0]));
202 memset(distances, 0, maxIdx * sizeof(distances[0]));
203 _kRingInternal(origin, k, out, distances, maxIdx, 0);
204 }
205}
206
207/**
208 * Internal helper function called recursively for kRingDistances.
209 *
210 * Adds the origin index to the output set (treating it as a hash set)
211 * and recurses to its neighbors, if needed.
212 *
213 * @param origin
214 * @param k Maximum distance to move from the origin.
215 * @param out Array treated as a hash set, elements being either H3Index or 0.
216 * @param distances Scratch area, with elements paralleling the out array.
217 * Elements indicate ijk distance from the origin index to the output index.
218 * @param maxIdx Size of out and scratch arrays (must be maxKringSize(k))
219 * @param curK Current distance from the origin.
220 */
221void _kRingInternal(H3Index origin, int k, H3Index* out, int* distances,
222 int maxIdx, int curK) {
223 if (origin == 0) return;
224
225 // Put origin in the output array. out is used as a hash set.
226 int off = origin % maxIdx;
227 while (out[off] != 0 && out[off] != origin) {
228 off = (off + 1) % maxIdx;
229 }
230
231 // We either got a free slot in the hash set or hit a duplicate
232 // We might need to process the duplicate anyways because we got
233 // here on a longer path before.
234 if (out[off] == origin && distances[off] <= curK) return;
235
236 out[off] = origin;
237 distances[off] = curK;
238
239 // Base case: reached an index k away from the origin.
240 if (curK >= k) return;
241
242 // Recurse to all neighbors in no particular order.
243 for (int i = 0; i < 6; i++) {
244 int rotations = 0;
245 _kRingInternal(h3NeighborRotations(origin, DIRECTIONS[i], &rotations),
246 k, out, distances, maxIdx, curK + 1);
247 }
248}
249
250/**
251 * Returns the hexagon index neighboring the origin, in the direction dir.
252 *
253 * Implementation note: The only reachable case where this returns 0 is if the
254 * origin is a pentagon and the translation is in the k direction. Thus,
255 * 0 can only be returned if origin is a pentagon.
256 *
257 * @param origin Origin index
258 * @param dir Direction to move in
259 * @param rotations Number of ccw rotations to perform to reorient the
260 * translation vector. Will be modified to the new number of
261 * rotations to perform (such as when crossing a face edge.)
262 * @return H3Index of the specified neighbor or 0 if deleted k-subsequence
263 * distortion is encountered.
264 */
265H3Index h3NeighborRotations(H3Index origin, Direction dir, int* rotations) {
266 H3Index out = origin;
267
268 for (int i = 0; i < *rotations; i++) {
269 dir = _rotate60ccw(dir);
270 }
271
272 int newRotations = 0;
273 int oldBaseCell = H3_GET_BASE_CELL(out);
274 Direction oldLeadingDigit = _h3LeadingNonZeroDigit(out);
275
276 // Adjust the indexing digits and, if needed, the base cell.
277 int r = H3_GET_RESOLUTION(out) - 1;
278 while (true) {
279 if (r == -1) {
280 H3_SET_BASE_CELL(out, baseCellNeighbors[oldBaseCell][dir]);
281 newRotations = baseCellNeighbor60CCWRots[oldBaseCell][dir];
282
283 if (H3_GET_BASE_CELL(out) == INVALID_BASE_CELL) {
284 // Adjust for the deleted k vertex at the base cell level.
285 // This edge actually borders a different neighbor.
286 H3_SET_BASE_CELL(out,
287 baseCellNeighbors[oldBaseCell][IK_AXES_DIGIT]);
288 newRotations =
289 baseCellNeighbor60CCWRots[oldBaseCell][IK_AXES_DIGIT];
290
291 // perform the adjustment for the k-subsequence we're skipping
292 // over.
293 out = _h3Rotate60ccw(out);
294 *rotations = *rotations + 1;
295 }
296
297 break;
298 } else {
299 Direction oldDigit = H3_GET_INDEX_DIGIT(out, r + 1);
300 Direction nextDir;
301 if (isResClassIII(r + 1)) {
302 H3_SET_INDEX_DIGIT(out, r + 1, NEW_DIGIT_II[oldDigit][dir]);
303 nextDir = NEW_ADJUSTMENT_II[oldDigit][dir];
304 } else {
305 H3_SET_INDEX_DIGIT(out, r + 1, NEW_DIGIT_III[oldDigit][dir]);
306 nextDir = NEW_ADJUSTMENT_III[oldDigit][dir];
307 }
308
309 if (nextDir != CENTER_DIGIT) {
310 dir = nextDir;
311 r--;
312 } else {
313 // No more adjustment to perform
314 break;
315 }
316 }
317 }
318
319 int newBaseCell = H3_GET_BASE_CELL(out);
320 if (_isBaseCellPentagon(newBaseCell)) {
321 int alreadyAdjustedKSubsequence = 0;
322
323 // force rotation out of missing k-axes sub-sequence
324 if (_h3LeadingNonZeroDigit(out) == K_AXES_DIGIT) {
325 if (oldBaseCell != newBaseCell) {
326 // in this case, we traversed into the deleted
327 // k subsequence of a pentagon base cell.
328 // We need to rotate out of that case depending
329 // on how we got here.
330 // check for a cw/ccw offset face; default is ccw
331
332 if (_baseCellIsCwOffset(
333 newBaseCell, baseCellData[oldBaseCell].homeFijk.face)) {
334 out = _h3Rotate60cw(out);
335 } else {
336 // See cwOffsetPent in testKRing.c for why this is
337 // unreachable.
338 out = _h3Rotate60ccw(out); // LCOV_EXCL_LINE
339 }
340 alreadyAdjustedKSubsequence = 1;
341 } else {
342 // In this case, we traversed into the deleted
343 // k subsequence from within the same pentagon
344 // base cell.
345 if (oldLeadingDigit == CENTER_DIGIT) {
346 // Undefined: the k direction is deleted from here
347 return H3_INVALID_INDEX;
348 } else if (oldLeadingDigit == JK_AXES_DIGIT) {
349 // Rotate out of the deleted k subsequence
350 // We also need an additional change to the direction we're
351 // moving in
352 out = _h3Rotate60ccw(out);
353 *rotations = *rotations + 1;
354 } else if (oldLeadingDigit == IK_AXES_DIGIT) {
355 // Rotate out of the deleted k subsequence
356 // We also need an additional change to the direction we're
357 // moving in
358 out = _h3Rotate60cw(out);
359 *rotations = *rotations + 5;
360 } else {
361 // Should never occur
362 return H3_INVALID_INDEX; // LCOV_EXCL_LINE
363 }
364 }
365 }
366
367 for (int i = 0; i < newRotations; i++) out = _h3RotatePent60ccw(out);
368
369 // Account for differing orientation of the base cells (this edge
370 // might not follow properties of some other edges.)
371 if (oldBaseCell != newBaseCell) {
372 if (_isBaseCellPolarPentagon(newBaseCell)) {
373 // 'polar' base cells behave differently because they have all
374 // i neighbors.
375 if (oldBaseCell != 118 && oldBaseCell != 8 &&
376 _h3LeadingNonZeroDigit(out) != JK_AXES_DIGIT) {
377 *rotations = *rotations + 1;
378 }
379 } else if (_h3LeadingNonZeroDigit(out) == IK_AXES_DIGIT &&
380 !alreadyAdjustedKSubsequence) {
381 // account for distortion introduced to the 5 neighbor by the
382 // deleted k subsequence.
383 *rotations = *rotations + 1;
384 }
385 }
386 } else {
387 for (int i = 0; i < newRotations; i++) out = _h3Rotate60ccw(out);
388 }
389
390 *rotations = (*rotations + newRotations) % 6;
391
392 return out;
393}
394
395/**
396 * hexRange produces indexes within k distance of the origin index.
397 * Output behavior is undefined when one of the indexes returned by this
398 * function is a pentagon or is in the pentagon distortion area.
399 *
400 * k-ring 0 is defined as the origin index, k-ring 1 is defined as k-ring 0 and
401 * all neighboring indexes, and so on.
402 *
403 * Output is placed in the provided array in order of increasing distance from
404 * the origin.
405 *
406 * @param origin Origin location.
407 * @param k k >= 0
408 * @param out Array which must be of size maxKringSize(k).
409 * @return 0 if no pentagon or pentagonal distortion area was encountered.
410 */
411int H3_EXPORT(hexRange)(H3Index origin, int k, H3Index* out) {
412 return H3_EXPORT(hexRangeDistances)(origin, k, out, 0);
413}
414
415/**
416 * hexRange produces indexes within k distance of the origin index.
417 * Output behavior is undefined when one of the indexes returned by this
418 * function is a pentagon or is in the pentagon distortion area.
419 *
420 * k-ring 0 is defined as the origin index, k-ring 1 is defined as k-ring 0 and
421 * all neighboring indexes, and so on.
422 *
423 * Output is placed in the provided array in order of increasing distance from
424 * the origin. The distances in hexagons is placed in the distances array at
425 * the same offset.
426 *
427 * @param origin Origin location.
428 * @param k k >= 0
429 * @param out Array which must be of size maxKringSize(k).
430 * @param distances Null or array which must be of size maxKringSize(k).
431 * @return 0 if no pentagon or pentagonal distortion area was encountered.
432 */
433int H3_EXPORT(hexRangeDistances)(H3Index origin, int k, H3Index* out,
434 int* distances) {
435 // Return codes:
436 // 1 Pentagon was encountered
437 // 2 Pentagon distortion (deleted k subsequence) was encountered
438 // Pentagon being encountered is not itself a problem; really the deleted
439 // k-subsequence is the problem, but for compatibility reasons we fail on
440 // the pentagon.
441
442 // k must be >= 0, so origin is always needed
443 int idx = 0;
444 out[idx] = origin;
445 if (distances) {
446 distances[idx] = 0;
447 }
448 idx++;
449
450 if (H3_EXPORT(h3IsPentagon)(origin)) {
451 // Pentagon was encountered; bail out as user doesn't want this.
452 return HEX_RANGE_PENTAGON;
453 }
454
455 // 0 < ring <= k, current ring
456 int ring = 1;
457 // 0 <= direction < 6, current side of the ring
458 int direction = 0;
459 // 0 <= i < ring, current position on the side of the ring
460 int i = 0;
461 // Number of 60 degree ccw rotations to perform on the direction (based on
462 // which faces have been crossed.)
463 int rotations = 0;
464
465 while (ring <= k) {
466 if (direction == 0 && i == 0) {
467 // Not putting in the output set as it will be done later, at
468 // the end of this ring.
469 origin =
470 h3NeighborRotations(origin, NEXT_RING_DIRECTION, &rotations);
471 if (origin == 0) {
472 // Should not be possible because `origin` would have to be a
473 // pentagon
474 return HEX_RANGE_K_SUBSEQUENCE; // LCOV_EXCL_LINE
475 }
476
477 if (H3_EXPORT(h3IsPentagon)(origin)) {
478 // Pentagon was encountered; bail out as user doesn't want this.
479 return HEX_RANGE_PENTAGON;
480 }
481 }
482
483 origin = h3NeighborRotations(origin, DIRECTIONS[direction], &rotations);
484 if (origin == 0) {
485 // Should not be possible because `origin` would have to be a
486 // pentagon
487 return HEX_RANGE_K_SUBSEQUENCE; // LCOV_EXCL_LINE
488 }
489 out[idx] = origin;
490 if (distances) {
491 distances[idx] = ring;
492 }
493 idx++;
494
495 i++;
496 // Check if end of this side of the k-ring
497 if (i == ring) {
498 i = 0;
499 direction++;
500 // Check if end of this ring.
501 if (direction == 6) {
502 direction = 0;
503 ring++;
504 }
505 }
506
507 if (H3_EXPORT(h3IsPentagon)(origin)) {
508 // Pentagon was encountered; bail out as user doesn't want this.
509 return HEX_RANGE_PENTAGON;
510 }
511 }
512 return HEX_RANGE_SUCCESS;
513}
514
515/**
516 * hexRanges takes an array of input hex IDs and a max k-ring and returns an
517 * array of hexagon IDs sorted first by the original hex IDs and then by the
518 * k-ring (0 to max), with no guaranteed sorting within each k-ring group.
519 *
520 * @param h3Set A pointer to an array of H3Indexes
521 * @param length The total number of H3Indexes in h3Set
522 * @param k The number of rings to generate
523 * @param out A pointer to the output memory to dump the new set of H3Indexes to
524 * The memory block should be equal to maxKringSize(k) * length
525 * @return 0 if no pentagon is encountered. Cannot trust output otherwise
526 */
527int H3_EXPORT(hexRanges)(H3Index* h3Set, int length, int k, H3Index* out) {
528 int success = 0;
529 H3Index* segment;
530 int segmentSize = H3_EXPORT(maxKringSize)(k);
531 for (int i = 0; i < length; i++) {
532 // Determine the appropriate segment of the output array to operate on
533 segment = out + i * segmentSize;
534 success = H3_EXPORT(hexRange)(h3Set[i], k, segment);
535 if (success != 0) return success;
536 }
537 return 0;
538}
539
540/**
541 * Returns the hollow hexagonal ring centered at origin with sides of length k.
542 *
543 * @param origin Origin location.
544 * @param k k >= 0
545 * @param out Array which must be of size 6 * k (or 1 if k == 0)
546 * @return 0 if no pentagonal distortion was encountered.
547 */
548int H3_EXPORT(hexRing)(H3Index origin, int k, H3Index* out) {
549 // Short-circuit on 'identity' ring
550 if (k == 0) {
551 out[0] = origin;
552 return 0;
553 }
554 int idx = 0;
555 // Number of 60 degree ccw rotations to perform on the direction (based on
556 // which faces have been crossed.)
557 int rotations = 0;
558 // Scratch structure for checking for pentagons
559 if (H3_EXPORT(h3IsPentagon)(origin)) {
560 // Pentagon was encountered; bail out as user doesn't want this.
561 return HEX_RANGE_PENTAGON;
562 }
563
564 for (int ring = 0; ring < k; ring++) {
565 origin = h3NeighborRotations(origin, NEXT_RING_DIRECTION, &rotations);
566 if (origin == 0) {
567 // Should not be possible because `origin` would have to be a
568 // pentagon
569 return HEX_RANGE_K_SUBSEQUENCE; // LCOV_EXCL_LINE
570 }
571
572 if (H3_EXPORT(h3IsPentagon)(origin)) {
573 return HEX_RANGE_PENTAGON;
574 }
575 }
576
577 H3Index lastIndex = origin;
578
579 out[idx] = origin;
580 idx++;
581
582 for (int direction = 0; direction < 6; direction++) {
583 for (int pos = 0; pos < k; pos++) {
584 origin =
585 h3NeighborRotations(origin, DIRECTIONS[direction], &rotations);
586 if (origin == 0) {
587 // Should not be possible because `origin` would have to be a
588 // pentagon
589 return HEX_RANGE_K_SUBSEQUENCE; // LCOV_EXCL_LINE
590 }
591
592 // Skip the very last index, it was already added. We do
593 // however need to traverse to it because of the pentagonal
594 // distortion check, below.
595 if (pos != k - 1 || direction != 5) {
596 out[idx] = origin;
597 idx++;
598
599 if (H3_EXPORT(h3IsPentagon)(origin)) {
600 return HEX_RANGE_PENTAGON;
601 }
602 }
603 }
604 }
605
606 // Check that this matches the expected lastIndex, if it doesn't,
607 // it indicates pentagonal distortion occurred and we should report
608 // failure.
609 if (lastIndex != origin) {
610 return HEX_RANGE_PENTAGON;
611 } else {
612 return HEX_RANGE_SUCCESS;
613 }
614}
615
616/**
617 * maxPolyfillSize returns the number of hexagons to allocate space for when
618 * performing a polyfill on the given GeoJSON-like data structure.
619 *
620 * Currently a laughably padded response, being a k-ring that wholly contains
621 * a bounding box of the GeoJSON, but still less wasted memory than initializing
622 * a Python application? ;)
623 *
624 * @param geoPolygon A GeoJSON-like data structure indicating the poly to fill
625 * @param res Hexagon resolution (0-15)
626 * @return number of hexagons to allocate for
627 */
628int H3_EXPORT(maxPolyfillSize)(const GeoPolygon* geoPolygon, int res) {
629 // Get the bounding box for the GeoJSON-like struct
630 BBox bbox;
631 bboxFromGeofence(&geoPolygon->geofence, &bbox);
632 int minK = bboxHexRadius(&bbox, res);
633
634 // The total number of hexagons to allocate can now be determined by
635 // the k-ring hex allocation helper function.
636 return H3_EXPORT(maxKringSize)(minK);
637}
638
639/**
640 * polyfill takes a given GeoJSON-like data structure and preallocated,
641 * zeroed memory, and fills it with the hexagons that are contained by
642 * the GeoJSON-like data structure.
643 *
644 * The current implementation is very primitive and slow, but correct,
645 * performing a point-in-poly operation on every hexagon in a k-ring defined
646 * around the given geofence.
647 *
648 * @param geoPolygon The geofence and holes defining the relevant area
649 * @param res The Hexagon resolution (0-15)
650 * @param out The slab of zeroed memory to write to. Assumed to be big enough.
651 */
652void H3_EXPORT(polyfill)(const GeoPolygon* geoPolygon, int res, H3Index* out) {
653 // One of the goals of the polyfill algorithm is that two adjacent polygons
654 // with zero overlap have zero overlapping hexagons. That the hexagons are
655 // uniquely assigned. There are a few approaches to take here, such as
656 // deciding based on which polygon has the greatest overlapping area of the
657 // hexagon, or the most number of contained points on the hexagon (using the
658 // center point as a tiebreaker).
659 //
660 // But if the polygons are convex, both of these more complex algorithms can
661 // be reduced down to checking whether or not the center of the hexagon is
662 // contained in the polygon, and so this is the approach that this polyfill
663 // algorithm will follow, as it's simpler, faster, and the error for concave
664 // polygons is still minimal (only affecting concave shapes on the order of
665 // magnitude of the hexagon size or smaller, not impacting larger concave
666 // shapes)
667 //
668 // This first part is identical to the maxPolyfillSize above.
669
670 // Get the bounding boxes for the polygon and any holes
671 BBox* bboxes = malloc((geoPolygon->numHoles + 1) * sizeof(BBox));
672 assert(bboxes != NULL);
673 bboxesFromGeoPolygon(geoPolygon, bboxes);
674 int minK = bboxHexRadius(&bboxes[0], res);
675 int numHexagons = H3_EXPORT(maxKringSize)(minK);
676
677 // Get the center hex
678 GeoCoord center;
679 bboxCenter(&bboxes[0], &center);
680 H3Index centerH3 = H3_EXPORT(geoToH3)(&center, res);
681
682 // From here on it works differently, first we get all potential
683 // hexagons inserted into the available memory
684 H3_EXPORT(kRing)(centerH3, minK, out);
685
686 // Next we iterate through each hexagon, and test its center point to see if
687 // it's contained in the GeoJSON-like struct
688 for (int i = 0; i < numHexagons; i++) {
689 // Skip records that are already zeroed
690 if (out[i] == 0) {
691 continue;
692 }
693 // Check if hexagon is inside of polygon
694 GeoCoord hexCenter;
695 H3_EXPORT(h3ToGeo)(out[i], &hexCenter);
696 hexCenter.lat = constrainLat(hexCenter.lat);
697 hexCenter.lon = constrainLng(hexCenter.lon);
698 // And remove from list if not
699 if (!pointInsidePolygon(geoPolygon, bboxes, &hexCenter)) {
700 out[i] = H3_INVALID_INDEX;
701 }
702 }
703 free(bboxes);
704}
705
706/**
707 * Internal: Create a vertex graph from a set of hexagons. It is the
708 * responsibility of the caller to call destroyVertexGraph on the populated
709 * graph, otherwise the memory in the graph nodes will not be freed.
710 * @private
711 * @param h3Set Set of hexagons
712 * @param numHexes Number of hexagons in the set
713 * @param graph Output graph
714 */
715void h3SetToVertexGraph(const H3Index* h3Set, const int numHexes,
716 VertexGraph* graph) {
717 GeoBoundary vertices;
718 GeoCoord* fromVtx;
719 GeoCoord* toVtx;
720 VertexNode* edge;
721 if (numHexes < 1) {
722 // We still need to init the graph, or calls to destroyVertexGraph will
723 // fail
724 initVertexGraph(graph, 0, 0);
725 return;
726 }
727 int res = H3_GET_RESOLUTION(h3Set[0]);
728 const int minBuckets = 6;
729 // TODO: Better way to calculate/guess?
730 int numBuckets = numHexes > minBuckets ? numHexes : minBuckets;
731 initVertexGraph(graph, numBuckets, res);
732 // Iterate through every hexagon
733 for (int i = 0; i < numHexes; i++) {
734 H3_EXPORT(h3ToGeoBoundary)(h3Set[i], &vertices);
735 // iterate through every edge
736 for (int j = 0; j < vertices.numVerts; j++) {
737 fromVtx = &vertices.verts[j];
738 toVtx = &vertices.verts[(j + 1) % vertices.numVerts];
739 // If we've seen this edge already, it will be reversed
740 edge = findNodeForEdge(graph, toVtx, fromVtx);
741 if (edge != NULL) {
742 // If we've seen it, drop it. No edge is shared by more than 2
743 // hexagons, so we'll never see it again.
744 removeVertexNode(graph, edge);
745 } else {
746 // Add a new node for this edge
747 addVertexNode(graph, fromVtx, toVtx);
748 }
749 }
750 }
751}
752
753/**
754 * Internal: Create a LinkedGeoPolygon from a vertex graph. It is the
755 * responsibility of the caller to call destroyLinkedPolygon on the populated
756 * linked geo structure, or the memory for that structure will not be freed.
757 * @private
758 * @param graph Input graph
759 * @param out Output polygon
760 */
761void _vertexGraphToLinkedGeo(VertexGraph* graph, LinkedGeoPolygon* out) {
762 *out = (LinkedGeoPolygon){0};
763 LinkedGeoLoop* loop;
764 VertexNode* edge;
765 GeoCoord nextVtx;
766 // Find the next unused entry point
767 while ((edge = firstVertexNode(graph)) != NULL) {
768 loop = addNewLinkedLoop(out);
769 // Walk the graph to get the outline
770 do {
771 addLinkedCoord(loop, &edge->from);
772 nextVtx = edge->to;
773 // Remove frees the node, so we can't use edge after this
774 removeVertexNode(graph, edge);
775 edge = findNodeForVertex(graph, &nextVtx);
776 } while (edge);
777 }
778}
779
780/**
781 * Create a LinkedGeoPolygon describing the outline(s) of a set of hexagons.
782 * Polygon outlines will follow GeoJSON MultiPolygon order: Each polygon will
783 * have one outer loop, which is first in the list, followed by any holes.
784 *
785 * It is the responsibility of the caller to call destroyLinkedPolygon on the
786 * populated linked geo structure, or the memory for that structure will
787 * not be freed.
788 *
789 * It is expected that all hexagons in the set have the same resolution and
790 * that the set contains no duplicates. Behavior is undefined if duplicates
791 * or multiple resolutions are present, and the algorithm may produce
792 * unexpected or invalid output.
793 *
794 * @param h3Set Set of hexagons
795 * @param numHexes Number of hexagons in set
796 * @param out Output polygon
797 */
798void H3_EXPORT(h3SetToLinkedGeo)(const H3Index* h3Set, const int numHexes,
799 LinkedGeoPolygon* out) {
800 VertexGraph graph;
801 h3SetToVertexGraph(h3Set, numHexes, &graph);
802 _vertexGraphToLinkedGeo(&graph, out);
803 // TODO: The return value, possibly indicating an error, is discarded here -
804 // we should use this when we update the API to return a value
805 normalizeMultiPolygon(out);
806 destroyVertexGraph(&graph);
807}
808