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 coordijk.c
17 * @brief Hex IJK coordinate systems functions including conversions to/from
18 * lat/lon.
19 */
20
21#include "coordijk.h"
22#include <math.h>
23#include <stdio.h>
24#include <stdlib.h>
25#include <string.h>
26#include "constants.h"
27#include "geoCoord.h"
28#include "mathExtensions.h"
29
30/**
31 * Sets an IJK coordinate to the specified component values.
32 *
33 * @param ijk The IJK coordinate to set.
34 * @param i The desired i component value.
35 * @param j The desired j component value.
36 * @param k The desired k component value.
37 */
38void _setIJK(CoordIJK* ijk, int i, int j, int k) {
39 ijk->i = i;
40 ijk->j = j;
41 ijk->k = k;
42}
43
44/**
45 * Determine the containing hex in ijk+ coordinates for a 2D cartesian
46 * coordinate vector (from DGGRID).
47 *
48 * @param v The 2D cartesian coordinate vector.
49 * @param h The ijk+ coordinates of the containing hex.
50 */
51void _hex2dToCoordIJK(const Vec2d* v, CoordIJK* h) {
52 double a1, a2;
53 double x1, x2;
54 int m1, m2;
55 double r1, r2;
56
57 // quantize into the ij system and then normalize
58 h->k = 0;
59
60 a1 = fabsl(v->x);
61 a2 = fabsl(v->y);
62
63 // first do a reverse conversion
64 x2 = a2 / M_SIN60;
65 x1 = a1 + x2 / 2.0L;
66
67 // check if we have the center of a hex
68 m1 = x1;
69 m2 = x2;
70
71 // otherwise round correctly
72 r1 = x1 - m1;
73 r2 = x2 - m2;
74
75 if (r1 < 0.5L) {
76 if (r1 < 1.0L / 3.0L) {
77 if (r2 < (1.0L + r1) / 2.0L) {
78 h->i = m1;
79 h->j = m2;
80 } else {
81 h->i = m1;
82 h->j = m2 + 1;
83 }
84 } else {
85 if (r2 < (1.0L - r1)) {
86 h->j = m2;
87 } else {
88 h->j = m2 + 1;
89 }
90
91 if ((1.0L - r1) <= r2 && r2 < (2.0 * r1)) {
92 h->i = m1 + 1;
93 } else {
94 h->i = m1;
95 }
96 }
97 } else {
98 if (r1 < 2.0L / 3.0L) {
99 if (r2 < (1.0L - r1)) {
100 h->j = m2;
101 } else {
102 h->j = m2 + 1;
103 }
104
105 if ((2.0L * r1 - 1.0L) < r2 && r2 < (1.0L - r1)) {
106 h->i = m1;
107 } else {
108 h->i = m1 + 1;
109 }
110 } else {
111 if (r2 < (r1 / 2.0L)) {
112 h->i = m1 + 1;
113 h->j = m2;
114 } else {
115 h->i = m1 + 1;
116 h->j = m2 + 1;
117 }
118 }
119 }
120
121 // now fold across the axes if necessary
122
123 if (v->x < 0.0L) {
124 if ((h->j % 2) == 0) // even
125 {
126 long long int axisi = h->j / 2;
127 long long int diff = h->i - axisi;
128 h->i = h->i - 2.0 * diff;
129 } else {
130 long long int axisi = (h->j + 1) / 2;
131 long long int diff = h->i - axisi;
132 h->i = h->i - (2.0 * diff + 1);
133 }
134 }
135
136 if (v->y < 0.0L) {
137 h->i = h->i - (2 * h->j + 1) / 2;
138 h->j = -1 * h->j;
139 }
140
141 _ijkNormalize(h);
142}
143
144/**
145 * Find the center point in 2D cartesian coordinates of a hex.
146 *
147 * @param h The ijk coordinates of the hex.
148 * @param v The 2D cartesian coordinates of the hex center point.
149 */
150void _ijkToHex2d(const CoordIJK* h, Vec2d* v) {
151 int i = h->i - h->k;
152 int j = h->j - h->k;
153
154 v->x = i - 0.5L * j;
155 v->y = j * M_SQRT3_2;
156}
157
158/**
159 * Returns whether or not two ijk coordinates contain exactly the same
160 * component values.
161 *
162 * @param c1 The first set of ijk coordinates.
163 * @param c2 The second set of ijk coordinates.
164 * @return 1 if the two addresses match, 0 if they do not.
165 */
166int _ijkMatches(const CoordIJK* c1, const CoordIJK* c2) {
167 return (c1->i == c2->i && c1->j == c2->j && c1->k == c2->k);
168}
169
170/**
171 * Add two ijk coordinates.
172 *
173 * @param h1 The first set of ijk coordinates.
174 * @param h2 The second set of ijk coordinates.
175 * @param sum The sum of the two sets of ijk coordinates.
176 */
177void _ijkAdd(const CoordIJK* h1, const CoordIJK* h2, CoordIJK* sum) {
178 sum->i = h1->i + h2->i;
179 sum->j = h1->j + h2->j;
180 sum->k = h1->k + h2->k;
181}
182
183/**
184 * Subtract two ijk coordinates.
185 *
186 * @param h1 The first set of ijk coordinates.
187 * @param h2 The second set of ijk coordinates.
188 * @param diff The difference of the two sets of ijk coordinates (h1 - h2).
189 */
190void _ijkSub(const CoordIJK* h1, const CoordIJK* h2, CoordIJK* diff) {
191 diff->i = h1->i - h2->i;
192 diff->j = h1->j - h2->j;
193 diff->k = h1->k - h2->k;
194}
195
196/**
197 * Uniformly scale ijk coordinates by a scalar. Works in place.
198 *
199 * @param c The ijk coordinates to scale.
200 * @param factor The scaling factor.
201 */
202void _ijkScale(CoordIJK* c, int factor) {
203 c->i *= factor;
204 c->j *= factor;
205 c->k *= factor;
206}
207
208/**
209 * Normalizes ijk coordinates by setting the components to the smallest possible
210 * values. Works in place.
211 *
212 * @param c The ijk coordinates to normalize.
213 */
214void _ijkNormalize(CoordIJK* c) {
215 // remove any negative values
216 if (c->i < 0) {
217 c->j -= c->i;
218 c->k -= c->i;
219 c->i = 0;
220 }
221
222 if (c->j < 0) {
223 c->i -= c->j;
224 c->k -= c->j;
225 c->j = 0;
226 }
227
228 if (c->k < 0) {
229 c->i -= c->k;
230 c->j -= c->k;
231 c->k = 0;
232 }
233
234 // remove the min value if needed
235 int min = c->i;
236 if (c->j < min) min = c->j;
237 if (c->k < min) min = c->k;
238 if (min > 0) {
239 c->i -= min;
240 c->j -= min;
241 c->k -= min;
242 }
243}
244
245/**
246 * Determines the H3 digit corresponding to a unit vector in ijk coordinates.
247 *
248 * @param ijk The ijk coordinates; must be a unit vector.
249 * @return The H3 digit (0-6) corresponding to the ijk unit vector, or
250 * INVALID_DIGIT on failure.
251 */
252Direction _unitIjkToDigit(const CoordIJK* ijk) {
253 CoordIJK c = *ijk;
254 _ijkNormalize(&c);
255
256 Direction digit = INVALID_DIGIT;
257 for (Direction i = CENTER_DIGIT; i < NUM_DIGITS; i++) {
258 if (_ijkMatches(&c, &UNIT_VECS[i])) {
259 digit = i;
260 break;
261 }
262 }
263
264 return digit;
265}
266
267/**
268 * Find the normalized ijk coordinates of the indexing parent of a cell in a
269 * counter-clockwise aperture 7 grid. Works in place.
270 *
271 * @param ijk The ijk coordinates.
272 */
273void _upAp7(CoordIJK* ijk) {
274 // convert to CoordIJ
275 int i = ijk->i - ijk->k;
276 int j = ijk->j - ijk->k;
277
278 ijk->i = (int)lroundl((3 * i - j) / 7.0L);
279 ijk->j = (int)lroundl((i + 2 * j) / 7.0L);
280 ijk->k = 0;
281 _ijkNormalize(ijk);
282}
283
284/**
285 * Find the normalized ijk coordinates of the indexing parent of a cell in a
286 * clockwise aperture 7 grid. Works in place.
287 *
288 * @param ijk The ijk coordinates.
289 */
290void _upAp7r(CoordIJK* ijk) {
291 // convert to CoordIJ
292 int i = ijk->i - ijk->k;
293 int j = ijk->j - ijk->k;
294
295 ijk->i = (int)lroundl((2 * i + j) / 7.0L);
296 ijk->j = (int)lroundl((3 * j - i) / 7.0L);
297 ijk->k = 0;
298 _ijkNormalize(ijk);
299}
300
301/**
302 * Find the normalized ijk coordinates of the hex centered on the indicated
303 * hex at the next finer aperture 7 counter-clockwise resolution. Works in
304 * place.
305 *
306 * @param ijk The ijk coordinates.
307 */
308void _downAp7(CoordIJK* ijk) {
309 // res r unit vectors in res r+1
310 CoordIJK iVec = {3, 0, 1};
311 CoordIJK jVec = {1, 3, 0};
312 CoordIJK kVec = {0, 1, 3};
313
314 _ijkScale(&iVec, ijk->i);
315 _ijkScale(&jVec, ijk->j);
316 _ijkScale(&kVec, ijk->k);
317
318 _ijkAdd(&iVec, &jVec, ijk);
319 _ijkAdd(ijk, &kVec, ijk);
320
321 _ijkNormalize(ijk);
322}
323
324/**
325 * Find the normalized ijk coordinates of the hex centered on the indicated
326 * hex at the next finer aperture 7 clockwise resolution. Works in place.
327 *
328 * @param ijk The ijk coordinates.
329 */
330void _downAp7r(CoordIJK* ijk) {
331 // res r unit vectors in res r+1
332 CoordIJK iVec = {3, 1, 0};
333 CoordIJK jVec = {0, 3, 1};
334 CoordIJK kVec = {1, 0, 3};
335
336 _ijkScale(&iVec, ijk->i);
337 _ijkScale(&jVec, ijk->j);
338 _ijkScale(&kVec, ijk->k);
339
340 _ijkAdd(&iVec, &jVec, ijk);
341 _ijkAdd(ijk, &kVec, ijk);
342
343 _ijkNormalize(ijk);
344}
345
346/**
347 * Find the normalized ijk coordinates of the hex in the specified digit
348 * direction from the specified ijk coordinates. Works in place.
349 *
350 * @param ijk The ijk coordinates.
351 * @param digit The digit direction from the original ijk coordinates.
352 */
353void _neighbor(CoordIJK* ijk, Direction digit) {
354 if (digit > CENTER_DIGIT && digit < NUM_DIGITS) {
355 _ijkAdd(ijk, &UNIT_VECS[digit], ijk);
356 _ijkNormalize(ijk);
357 }
358}
359
360/**
361 * Rotates ijk coordinates 60 degrees counter-clockwise. Works in place.
362 *
363 * @param ijk The ijk coordinates.
364 */
365void _ijkRotate60ccw(CoordIJK* ijk) {
366 // unit vector rotations
367 CoordIJK iVec = {1, 1, 0};
368 CoordIJK jVec = {0, 1, 1};
369 CoordIJK kVec = {1, 0, 1};
370
371 _ijkScale(&iVec, ijk->i);
372 _ijkScale(&jVec, ijk->j);
373 _ijkScale(&kVec, ijk->k);
374
375 _ijkAdd(&iVec, &jVec, ijk);
376 _ijkAdd(ijk, &kVec, ijk);
377
378 _ijkNormalize(ijk);
379}
380
381/**
382 * Rotates ijk coordinates 60 degrees clockwise. Works in place.
383 *
384 * @param ijk The ijk coordinates.
385 */
386void _ijkRotate60cw(CoordIJK* ijk) {
387 // unit vector rotations
388 CoordIJK iVec = {1, 0, 1};
389 CoordIJK jVec = {1, 1, 0};
390 CoordIJK kVec = {0, 1, 1};
391
392 _ijkScale(&iVec, ijk->i);
393 _ijkScale(&jVec, ijk->j);
394 _ijkScale(&kVec, ijk->k);
395
396 _ijkAdd(&iVec, &jVec, ijk);
397 _ijkAdd(ijk, &kVec, ijk);
398
399 _ijkNormalize(ijk);
400}
401
402/**
403 * Rotates indexing digit 60 degrees counter-clockwise. Returns result.
404 *
405 * @param digit Indexing digit (between 1 and 6 inclusive)
406 */
407Direction _rotate60ccw(Direction digit) {
408 switch (digit) {
409 case K_AXES_DIGIT:
410 return IK_AXES_DIGIT;
411 case IK_AXES_DIGIT:
412 return I_AXES_DIGIT;
413 case I_AXES_DIGIT:
414 return IJ_AXES_DIGIT;
415 case IJ_AXES_DIGIT:
416 return J_AXES_DIGIT;
417 case J_AXES_DIGIT:
418 return JK_AXES_DIGIT;
419 case JK_AXES_DIGIT:
420 return K_AXES_DIGIT;
421 default:
422 return digit;
423 }
424}
425
426/**
427 * Rotates indexing digit 60 degrees clockwise. Returns result.
428 *
429 * @param digit Indexing digit (between 1 and 6 inclusive)
430 */
431Direction _rotate60cw(Direction digit) {
432 switch (digit) {
433 case K_AXES_DIGIT:
434 return JK_AXES_DIGIT;
435 case JK_AXES_DIGIT:
436 return J_AXES_DIGIT;
437 case J_AXES_DIGIT:
438 return IJ_AXES_DIGIT;
439 case IJ_AXES_DIGIT:
440 return I_AXES_DIGIT;
441 case I_AXES_DIGIT:
442 return IK_AXES_DIGIT;
443 case IK_AXES_DIGIT:
444 return K_AXES_DIGIT;
445 default:
446 return digit;
447 }
448}
449
450/**
451 * Find the normalized ijk coordinates of the hex centered on the indicated
452 * hex at the next finer aperture 3 counter-clockwise resolution. Works in
453 * place.
454 *
455 * @param ijk The ijk coordinates.
456 */
457void _downAp3(CoordIJK* ijk) {
458 // res r unit vectors in res r+1
459 CoordIJK iVec = {2, 0, 1};
460 CoordIJK jVec = {1, 2, 0};
461 CoordIJK kVec = {0, 1, 2};
462
463 _ijkScale(&iVec, ijk->i);
464 _ijkScale(&jVec, ijk->j);
465 _ijkScale(&kVec, ijk->k);
466
467 _ijkAdd(&iVec, &jVec, ijk);
468 _ijkAdd(ijk, &kVec, ijk);
469
470 _ijkNormalize(ijk);
471}
472
473/**
474 * Find the normalized ijk coordinates of the hex centered on the indicated
475 * hex at the next finer aperture 3 clockwise resolution. Works in place.
476 *
477 * @param ijk The ijk coordinates.
478 */
479void _downAp3r(CoordIJK* ijk) {
480 // res r unit vectors in res r+1
481 CoordIJK iVec = {2, 1, 0};
482 CoordIJK jVec = {0, 2, 1};
483 CoordIJK kVec = {1, 0, 2};
484
485 _ijkScale(&iVec, ijk->i);
486 _ijkScale(&jVec, ijk->j);
487 _ijkScale(&kVec, ijk->k);
488
489 _ijkAdd(&iVec, &jVec, ijk);
490 _ijkAdd(ijk, &kVec, ijk);
491
492 _ijkNormalize(ijk);
493}
494
495/**
496 * Finds the distance between the two coordinates. Returns result.
497 *
498 * @param c1 The first set of ijk coordinates.
499 * @param c2 The second set of ijk coordinates.
500 */
501int ijkDistance(const CoordIJK* c1, const CoordIJK* c2) {
502 CoordIJK diff;
503 _ijkSub(c1, c2, &diff);
504 _ijkNormalize(&diff);
505 CoordIJK absDiff = {abs(diff.i), abs(diff.j), abs(diff.k)};
506 return MAX(absDiff.i, MAX(absDiff.j, absDiff.k));
507}
508
509/**
510 * Transforms coordinates from the IJK+ coordinate system to the IJ coordinate
511 * system.
512 *
513 * @param ijk The input IJK+ coordinates
514 * @param ij The output IJ coordinates
515 */
516void ijkToIj(const CoordIJK* ijk, CoordIJ* ij) {
517 ij->i = ijk->i - ijk->k;
518 ij->j = ijk->j - ijk->k;
519}
520
521/**
522 * Transforms coordinates from the IJ coordinate system to the IJK+ coordinate
523 * system.
524 *
525 * @param ij The input IJ coordinates
526 * @param ijk The output IJK+ coordinates
527 */
528void ijToIjk(const CoordIJ* ij, CoordIJK* ijk) {
529 ijk->i = ij->i;
530 ijk->j = ij->j;
531 ijk->k = 0;
532
533 _ijkNormalize(ijk);
534}
535
536/**
537 * Convert IJK coordinates to cube coordinates, in place
538 * @param ijk Coordinate to convert
539 */
540void ijkToCube(CoordIJK* ijk) {
541 ijk->i = -ijk->i + ijk->k;
542 ijk->j = ijk->j - ijk->k;
543 ijk->k = -ijk->i - ijk->j;
544}
545
546/**
547 * Convert cube coordinates to IJK coordinates, in place
548 * @param ijk Coordinate to convert
549 */
550void cubeToIjk(CoordIJK* ijk) {
551 ijk->i = -ijk->i;
552 ijk->k = 0;
553 _ijkNormalize(ijk);
554}