/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 | } |