1/*
2 * Copyright 2016-2017 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 faceijk.c
17 * @brief Functions for working with icosahedral face-centered hex IJK
18 * coordinate systems.
19 */
20
21#include "faceijk.h"
22#include <assert.h>
23#include <math.h>
24#include <stdio.h>
25#include <stdlib.h>
26#include <string.h>
27#include "constants.h"
28#include "coordijk.h"
29#include "geoCoord.h"
30#include "h3Index.h"
31#include "vec3d.h"
32
33/** square root of 7 */
34#define M_SQRT7 2.6457513110645905905016157536392604257102L
35
36/** @brief icosahedron face centers in lat/lon radians */
37const GeoCoord faceCenterGeo[NUM_ICOSA_FACES] = {
38 {0.803582649718989942, 1.248397419617396099}, // face 0
39 {1.307747883455638156, 2.536945009877921159}, // face 1
40 {1.054751253523952054, -1.347517358900396623}, // face 2
41 {0.600191595538186799, -0.450603909469755746}, // face 3
42 {0.491715428198773866, 0.401988202911306943}, // face 4
43 {0.172745327415618701, 1.678146885280433686}, // face 5
44 {0.605929321571350690, 2.953923329812411617}, // face 6
45 {0.427370518328979641, -1.888876200336285401}, // face 7
46 {-0.079066118549212831, -0.733429513380867741}, // face 8
47 {-0.230961644455383637, 0.506495587332349035}, // face 9
48 {0.079066118549212831, 2.408163140208925497}, // face 10
49 {0.230961644455383637, -2.635097066257444203}, // face 11
50 {-0.172745327415618701, -1.463445768309359553}, // face 12
51 {-0.605929321571350690, -0.187669323777381622}, // face 13
52 {-0.427370518328979641, 1.252716453253507838}, // face 14
53 {-0.600191595538186799, 2.690988744120037492}, // face 15
54 {-0.491715428198773866, -2.739604450678486295}, // face 16
55 {-0.803582649718989942, -1.893195233972397139}, // face 17
56 {-1.307747883455638156, -0.604647643711872080}, // face 18
57 {-1.054751253523952054, 1.794075294689396615}, // face 19
58};
59
60/** @brief icosahedron face centers in x/y/z on the unit sphere */
61static const Vec3d faceCenterPoint[NUM_ICOSA_FACES] = {
62 {0.2199307791404606, 0.6583691780274996, 0.7198475378926182}, // face 0
63 {-0.2139234834501421, 0.1478171829550703, 0.9656017935214205}, // face 1
64 {0.1092625278784797, -0.4811951572873210, 0.8697775121287253}, // face 2
65 {0.7428567301586791, -0.3593941678278028, 0.5648005936517033}, // face 3
66 {0.8112534709140969, 0.3448953237639384, 0.4721387736413930}, // face 4
67 {-0.1055498149613921, 0.9794457296411413, 0.1718874610009365}, // face 5
68 {-0.8075407579970092, 0.1533552485898818, 0.5695261994882688}, // face 6
69 {-0.2846148069787907, -0.8644080972654206, 0.4144792552473539}, // face 7
70 {0.7405621473854482, -0.6673299564565524, -0.0789837646326737}, // face 8
71 {0.8512303986474293, 0.4722343788582681, -0.2289137388687808}, // face 9
72 {-0.7405621473854481, 0.6673299564565524, 0.0789837646326737}, // face 10
73 {-0.8512303986474292, -0.4722343788582682, 0.2289137388687808}, // face 11
74 {0.1055498149613919, -0.9794457296411413, -0.1718874610009365}, // face 12
75 {0.8075407579970092, -0.1533552485898819, -0.5695261994882688}, // face 13
76 {0.2846148069787908, 0.8644080972654204, -0.4144792552473539}, // face 14
77 {-0.7428567301586791, 0.3593941678278027, -0.5648005936517033}, // face 15
78 {-0.8112534709140971, -0.3448953237639382, -0.4721387736413930}, // face 16
79 {-0.2199307791404607, -0.6583691780274996, -0.7198475378926182}, // face 17
80 {0.2139234834501420, -0.1478171829550704, -0.9656017935214205}, // face 18
81 {-0.1092625278784796, 0.4811951572873210, -0.8697775121287253}, // face 19
82};
83
84/** @brief icosahedron face ijk axes as azimuth in radians from face center to
85 * vertex 0/1/2 respectively
86 */
87static const double faceAxesAzRadsCII[NUM_ICOSA_FACES][3] = {
88 {5.619958268523939882, 3.525563166130744542,
89 1.431168063737548730}, // face 0
90 {5.760339081714187279, 3.665943979320991689,
91 1.571548876927796127}, // face 1
92 {0.780213654393430055, 4.969003859179821079,
93 2.874608756786625655}, // face 2
94 {0.430469363979999913, 4.619259568766391033,
95 2.524864466373195467}, // face 3
96 {6.130269123335111400, 4.035874020941915804,
97 1.941478918548720291}, // face 4
98 {2.692877706530642877, 0.598482604137447119,
99 4.787272808923838195}, // face 5
100 {2.982963003477243874, 0.888567901084048369,
101 5.077358105870439581}, // face 6
102 {3.532912002790141181, 1.438516900396945656,
103 5.627307105183336758}, // face 7
104 {3.494305004259568154, 1.399909901866372864,
105 5.588700106652763840}, // face 8
106 {3.003214169499538391, 0.908819067106342928,
107 5.097609271892733906}, // face 9
108 {5.930472956509811562, 3.836077854116615875,
109 1.741682751723420374}, // face 10
110 {0.138378484090254847, 4.327168688876645809,
111 2.232773586483450311}, // face 11
112 {0.448714947059150361, 4.637505151845541521,
113 2.543110049452346120}, // face 12
114 {0.158629650112549365, 4.347419854898940135,
115 2.253024752505744869}, // face 13
116 {5.891865957979238535, 3.797470855586042958,
117 1.703075753192847583}, // face 14
118 {2.711123289609793325, 0.616728187216597771,
119 4.805518392002988683}, // face 15
120 {3.294508837434268316, 1.200113735041072948,
121 5.388903939827463911}, // face 16
122 {3.804819692245439833, 1.710424589852244509,
123 5.899214794638635174}, // face 17
124 {3.664438879055192436, 1.570043776661997111,
125 5.758833981448388027}, // face 18
126 {2.361378999196363184, 0.266983896803167583,
127 4.455774101589558636}, // face 19
128};
129
130/** @brief Definition of which faces neighbor each other. */
131static const FaceOrientIJK faceNeighbors[NUM_ICOSA_FACES][4] = {
132 {
133 // face 0
134 {0, {0, 0, 0}, 0}, // central face
135 {4, {2, 0, 2}, 1}, // ij quadrant
136 {1, {2, 2, 0}, 5}, // ki quadrant
137 {5, {0, 2, 2}, 3} // jk quadrant
138 },
139 {
140 // face 1
141 {1, {0, 0, 0}, 0}, // central face
142 {0, {2, 0, 2}, 1}, // ij quadrant
143 {2, {2, 2, 0}, 5}, // ki quadrant
144 {6, {0, 2, 2}, 3} // jk quadrant
145 },
146 {
147 // face 2
148 {2, {0, 0, 0}, 0}, // central face
149 {1, {2, 0, 2}, 1}, // ij quadrant
150 {3, {2, 2, 0}, 5}, // ki quadrant
151 {7, {0, 2, 2}, 3} // jk quadrant
152 },
153 {
154 // face 3
155 {3, {0, 0, 0}, 0}, // central face
156 {2, {2, 0, 2}, 1}, // ij quadrant
157 {4, {2, 2, 0}, 5}, // ki quadrant
158 {8, {0, 2, 2}, 3} // jk quadrant
159 },
160 {
161 // face 4
162 {4, {0, 0, 0}, 0}, // central face
163 {3, {2, 0, 2}, 1}, // ij quadrant
164 {0, {2, 2, 0}, 5}, // ki quadrant
165 {9, {0, 2, 2}, 3} // jk quadrant
166 },
167 {
168 // face 5
169 {5, {0, 0, 0}, 0}, // central face
170 {10, {2, 2, 0}, 3}, // ij quadrant
171 {14, {2, 0, 2}, 3}, // ki quadrant
172 {0, {0, 2, 2}, 3} // jk quadrant
173 },
174 {
175 // face 6
176 {6, {0, 0, 0}, 0}, // central face
177 {11, {2, 2, 0}, 3}, // ij quadrant
178 {10, {2, 0, 2}, 3}, // ki quadrant
179 {1, {0, 2, 2}, 3} // jk quadrant
180 },
181 {
182 // face 7
183 {7, {0, 0, 0}, 0}, // central face
184 {12, {2, 2, 0}, 3}, // ij quadrant
185 {11, {2, 0, 2}, 3}, // ki quadrant
186 {2, {0, 2, 2}, 3} // jk quadrant
187 },
188 {
189 // face 8
190 {8, {0, 0, 0}, 0}, // central face
191 {13, {2, 2, 0}, 3}, // ij quadrant
192 {12, {2, 0, 2}, 3}, // ki quadrant
193 {3, {0, 2, 2}, 3} // jk quadrant
194 },
195 {
196 // face 9
197 {9, {0, 0, 0}, 0}, // central face
198 {14, {2, 2, 0}, 3}, // ij quadrant
199 {13, {2, 0, 2}, 3}, // ki quadrant
200 {4, {0, 2, 2}, 3} // jk quadrant
201 },
202 {
203 // face 10
204 {10, {0, 0, 0}, 0}, // central face
205 {5, {2, 2, 0}, 3}, // ij quadrant
206 {6, {2, 0, 2}, 3}, // ki quadrant
207 {15, {0, 2, 2}, 3} // jk quadrant
208 },
209 {
210 // face 11
211 {11, {0, 0, 0}, 0}, // central face
212 {6, {2, 2, 0}, 3}, // ij quadrant
213 {7, {2, 0, 2}, 3}, // ki quadrant
214 {16, {0, 2, 2}, 3} // jk quadrant
215 },
216 {
217 // face 12
218 {12, {0, 0, 0}, 0}, // central face
219 {7, {2, 2, 0}, 3}, // ij quadrant
220 {8, {2, 0, 2}, 3}, // ki quadrant
221 {17, {0, 2, 2}, 3} // jk quadrant
222 },
223 {
224 // face 13
225 {13, {0, 0, 0}, 0}, // central face
226 {8, {2, 2, 0}, 3}, // ij quadrant
227 {9, {2, 0, 2}, 3}, // ki quadrant
228 {18, {0, 2, 2}, 3} // jk quadrant
229 },
230 {
231 // face 14
232 {14, {0, 0, 0}, 0}, // central face
233 {9, {2, 2, 0}, 3}, // ij quadrant
234 {5, {2, 0, 2}, 3}, // ki quadrant
235 {19, {0, 2, 2}, 3} // jk quadrant
236 },
237 {
238 // face 15
239 {15, {0, 0, 0}, 0}, // central face
240 {16, {2, 0, 2}, 1}, // ij quadrant
241 {19, {2, 2, 0}, 5}, // ki quadrant
242 {10, {0, 2, 2}, 3} // jk quadrant
243 },
244 {
245 // face 16
246 {16, {0, 0, 0}, 0}, // central face
247 {17, {2, 0, 2}, 1}, // ij quadrant
248 {15, {2, 2, 0}, 5}, // ki quadrant
249 {11, {0, 2, 2}, 3} // jk quadrant
250 },
251 {
252 // face 17
253 {17, {0, 0, 0}, 0}, // central face
254 {18, {2, 0, 2}, 1}, // ij quadrant
255 {16, {2, 2, 0}, 5}, // ki quadrant
256 {12, {0, 2, 2}, 3} // jk quadrant
257 },
258 {
259 // face 18
260 {18, {0, 0, 0}, 0}, // central face
261 {19, {2, 0, 2}, 1}, // ij quadrant
262 {17, {2, 2, 0}, 5}, // ki quadrant
263 {13, {0, 2, 2}, 3} // jk quadrant
264 },
265 {
266 // face 19
267 {19, {0, 0, 0}, 0}, // central face
268 {15, {2, 0, 2}, 1}, // ij quadrant
269 {18, {2, 2, 0}, 5}, // ki quadrant
270 {14, {0, 2, 2}, 3} // jk quadrant
271 }};
272
273/** @brief direction from the origin face to the destination face, relative to
274 * the origin face's coordinate system, or -1 if not adjacent.
275 */
276static const int adjacentFaceDir[NUM_ICOSA_FACES][NUM_ICOSA_FACES] = {
277 {0, KI, -1, -1, IJ, JK, -1, -1, -1, -1,
278 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 0
279 {IJ, 0, KI, -1, -1, -1, JK, -1, -1, -1,
280 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 1
281 {-1, IJ, 0, KI, -1, -1, -1, JK, -1, -1,
282 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 2
283 {-1, -1, IJ, 0, KI, -1, -1, -1, JK, -1,
284 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 3
285 {KI, -1, -1, IJ, 0, -1, -1, -1, -1, JK,
286 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 4
287 {JK, -1, -1, -1, -1, 0, -1, -1, -1, -1,
288 IJ, -1, -1, -1, KI, -1, -1, -1, -1, -1}, // face 5
289 {-1, JK, -1, -1, -1, -1, 0, -1, -1, -1,
290 KI, IJ, -1, -1, -1, -1, -1, -1, -1, -1}, // face 6
291 {-1, -1, JK, -1, -1, -1, -1, 0, -1, -1,
292 -1, KI, IJ, -1, -1, -1, -1, -1, -1, -1}, // face 7
293 {-1, -1, -1, JK, -1, -1, -1, -1, 0, -1,
294 -1, -1, KI, IJ, -1, -1, -1, -1, -1, -1}, // face 8
295 {-1, -1, -1, -1, JK, -1, -1, -1, -1, 0,
296 -1, -1, -1, KI, IJ, -1, -1, -1, -1, -1}, // face 9
297 {-1, -1, -1, -1, -1, IJ, KI, -1, -1, -1,
298 0, -1, -1, -1, -1, JK, -1, -1, -1, -1}, // face 10
299 {-1, -1, -1, -1, -1, -1, IJ, KI, -1, -1,
300 -1, 0, -1, -1, -1, -1, JK, -1, -1, -1}, // face 11
301 {-1, -1, -1, -1, -1, -1, -1, IJ, KI, -1,
302 -1, -1, 0, -1, -1, -1, -1, JK, -1, -1}, // face 12
303 {-1, -1, -1, -1, -1, -1, -1, -1, IJ, KI,
304 -1, -1, -1, 0, -1, -1, -1, -1, JK, -1}, // face 13
305 {-1, -1, -1, -1, -1, KI, -1, -1, -1, IJ,
306 -1, -1, -1, -1, 0, -1, -1, -1, -1, JK}, // face 14
307 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
308 JK, -1, -1, -1, -1, 0, IJ, -1, -1, KI}, // face 15
309 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
310 -1, JK, -1, -1, -1, KI, 0, IJ, -1, -1}, // face 16
311 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
312 -1, -1, JK, -1, -1, -1, KI, 0, IJ, -1}, // face 17
313 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
314 -1, -1, -1, JK, -1, -1, -1, KI, 0, IJ}, // face 18
315 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
316 -1, -1, -1, -1, JK, IJ, -1, -1, KI, 0} // face 19
317};
318
319/** @brief overage distance table */
320static const int maxDimByCIIres[] = {
321 2, // res 0
322 -1, // res 1
323 14, // res 2
324 -1, // res 3
325 98, // res 4
326 -1, // res 5
327 686, // res 6
328 -1, // res 7
329 4802, // res 8
330 -1, // res 9
331 33614, // res 10
332 -1, // res 11
333 235298, // res 12
334 -1, // res 13
335 1647086, // res 14
336 -1, // res 15
337 11529602 // res 16
338};
339
340/** @brief unit scale distance table */
341static const int unitScaleByCIIres[] = {
342 1, // res 0
343 -1, // res 1
344 7, // res 2
345 -1, // res 3
346 49, // res 4
347 -1, // res 5
348 343, // res 6
349 -1, // res 7
350 2401, // res 8
351 -1, // res 9
352 16807, // res 10
353 -1, // res 11
354 117649, // res 12
355 -1, // res 13
356 823543, // res 14
357 -1, // res 15
358 5764801 // res 16
359};
360
361/**
362 * Encodes a coordinate on the sphere to the FaceIJK address of the containing
363 * cell at the specified resolution.
364 *
365 * @param g The spherical coordinates to encode.
366 * @param res The desired H3 resolution for the encoding.
367 * @param h The FaceIJK address of the containing cell at resolution res.
368 */
369void _geoToFaceIjk(const GeoCoord* g, int res, FaceIJK* h) {
370 // first convert to hex2d
371 Vec2d v;
372 _geoToHex2d(g, res, &h->face, &v);
373
374 // then convert to ijk+
375 _hex2dToCoordIJK(&v, &h->coord);
376}
377
378/**
379 * Encodes a coordinate on the sphere to the corresponding icosahedral face and
380 * containing 2D hex coordinates relative to that face center.
381 *
382 * @param g The spherical coordinates to encode.
383 * @param res The desired H3 resolution for the encoding.
384 * @param face The icosahedral face containing the spherical coordinates.
385 * @param v The 2D hex coordinates of the cell containing the point.
386 */
387void _geoToHex2d(const GeoCoord* g, int res, int* face, Vec2d* v) {
388 Vec3d v3d;
389 _geoToVec3d(g, &v3d);
390
391 // determine the icosahedron face
392 *face = 0;
393 double sqd = _pointSquareDist(&faceCenterPoint[0], &v3d);
394 for (int f = 1; f < NUM_ICOSA_FACES; f++) {
395 double sqdT = _pointSquareDist(&faceCenterPoint[f], &v3d);
396 if (sqdT < sqd) {
397 *face = f;
398 sqd = sqdT;
399 }
400 }
401
402 // cos(r) = 1 - 2 * sin^2(r/2) = 1 - 2 * (sqd / 4) = 1 - sqd/2
403 double r = acos(1 - sqd / 2);
404
405 if (r < EPSILON) {
406 v->x = v->y = 0.0L;
407 return;
408 }
409
410 // now have face and r, now find CCW theta from CII i-axis
411 double theta =
412 _posAngleRads(faceAxesAzRadsCII[*face][0] -
413 _posAngleRads(_geoAzimuthRads(&faceCenterGeo[*face], g)));
414
415 // adjust theta for Class III (odd resolutions)
416 if (isResClassIII(res)) theta = _posAngleRads(theta - M_AP7_ROT_RADS);
417
418 // perform gnomonic scaling of r
419 r = tan(r);
420
421 // scale for current resolution length u
422 r /= RES0_U_GNOMONIC;
423 for (int i = 0; i < res; i++) r *= M_SQRT7;
424
425 // we now have (r, theta) in hex2d with theta ccw from x-axes
426
427 // convert to local x,y
428 v->x = r * cos(theta);
429 v->y = r * sin(theta);
430}
431
432/**
433 * Determines the center point in spherical coordinates of a cell given by 2D
434 * hex coordinates on a particular icosahedral face.
435 *
436 * @param v The 2D hex coordinates of the cell.
437 * @param face The icosahedral face upon which the 2D hex coordinate system is
438 * centered.
439 * @param res The H3 resolution of the cell.
440 * @param substrate Indicates whether or not this grid is actually a substrate
441 * grid relative to the specified resolution.
442 * @param g The spherical coordinates of the cell center point.
443 */
444void _hex2dToGeo(const Vec2d* v, int face, int res, int substrate,
445 GeoCoord* g) {
446 // calculate (r, theta) in hex2d
447 double r = _v2dMag(v);
448
449 if (r < EPSILON) {
450 *g = faceCenterGeo[face];
451 return;
452 }
453
454 double theta = atan2(v->y, v->x);
455
456 // scale for current resolution length u
457 for (int i = 0; i < res; i++) r /= M_SQRT7;
458
459 // scale accordingly if this is a substrate grid
460 if (substrate) {
461 r /= 3.0;
462 if (isResClassIII(res)) r /= M_SQRT7;
463 }
464
465 r *= RES0_U_GNOMONIC;
466
467 // perform inverse gnomonic scaling of r
468 r = atan(r);
469
470 // adjust theta for Class III
471 // if a substrate grid, then it's already been adjusted for Class III
472 if (!substrate && isResClassIII(res))
473 theta = _posAngleRads(theta + M_AP7_ROT_RADS);
474
475 // find theta as an azimuth
476 theta = _posAngleRads(faceAxesAzRadsCII[face][0] - theta);
477
478 // now find the point at (r,theta) from the face center
479 _geoAzDistanceRads(&faceCenterGeo[face], theta, r, g);
480}
481
482/**
483 * Determines the center point in spherical coordinates of a cell given by
484 * a FaceIJK address at a specified resolution.
485 *
486 * @param h The FaceIJK address of the cell.
487 * @param res The H3 resolution of the cell.
488 * @param g The spherical coordinates of the cell center point.
489 */
490void _faceIjkToGeo(const FaceIJK* h, int res, GeoCoord* g) {
491 Vec2d v;
492 _ijkToHex2d(&h->coord, &v);
493 _hex2dToGeo(&v, h->face, res, 0, g);
494}
495
496/**
497 * Generates the cell boundary in spherical coordinates for a pentagonal cell
498 * given by a FaceIJK address at a specified resolution.
499 *
500 * @param h The FaceIJK address of the pentagonal cell.
501 * @param res The H3 resolution of the cell.
502 * @param g The spherical coordinates of the cell boundary.
503 */
504void _faceIjkPentToGeoBoundary(const FaceIJK* h, int res, GeoBoundary* g) {
505 // the vertexes of an origin-centered pentagon in a Class II resolution on a
506 // substrate grid with aperture sequence 33r. The aperture 3 gets us the
507 // vertices, and the 3r gets us back to Class II.
508 // vertices listed ccw from the i-axes
509 CoordIJK vertsCII[NUM_PENT_VERTS] = {
510 {2, 1, 0}, // 0
511 {1, 2, 0}, // 1
512 {0, 2, 1}, // 2
513 {0, 1, 2}, // 3
514 {1, 0, 2}, // 4
515 };
516
517 // the vertexes of an origin-centered pentagon in a Class III resolution on
518 // a substrate grid with aperture sequence 33r7r. The aperture 3 gets us the
519 // vertices, and the 3r7r gets us to Class II. vertices listed ccw from the
520 // i-axes
521 CoordIJK vertsCIII[NUM_PENT_VERTS] = {
522 {5, 4, 0}, // 0
523 {1, 5, 0}, // 1
524 {0, 5, 4}, // 2
525 {0, 1, 5}, // 3
526 {4, 0, 5}, // 4
527 };
528
529 // get the correct set of substrate vertices for this resolution
530 CoordIJK* verts;
531 if (isResClassIII(res))
532 verts = vertsCIII;
533 else
534 verts = vertsCII;
535
536 // adjust the center point to be in an aperture 33r substrate grid
537 // these should be composed for speed
538 FaceIJK centerIJK = *h;
539 _downAp3(&centerIJK.coord);
540 _downAp3r(&centerIJK.coord);
541
542 // if res is Class III we need to add a cw aperture 7 to get to
543 // icosahedral Class II
544 int adjRes = res;
545 if (isResClassIII(res)) {
546 _downAp7r(&centerIJK.coord);
547 adjRes++;
548 }
549
550 // The center point is now in the same substrate grid as the origin
551 // cell vertices. Add the center point substate coordinates
552 // to each vertex to translate the vertices to that cell.
553 FaceIJK fijkVerts[NUM_PENT_VERTS];
554 for (int v = 0; v < NUM_PENT_VERTS; v++) {
555 fijkVerts[v].face = centerIJK.face;
556 _ijkAdd(&centerIJK.coord, &verts[v], &fijkVerts[v].coord);
557 _ijkNormalize(&fijkVerts[v].coord);
558 }
559
560 // convert each vertex to lat/lon
561 // adjust the face of each vertex as appropriate and introduce
562 // edge-crossing vertices as needed
563 g->numVerts = 0;
564 FaceIJK lastFijk;
565 for (int vert = 0; vert < NUM_PENT_VERTS + 1; vert++) {
566 int v = vert % NUM_PENT_VERTS;
567
568 FaceIJK fijk = fijkVerts[v];
569
570 int pentLeading4 = 0;
571 int overage = _adjustOverageClassII(&fijk, adjRes, pentLeading4, 1);
572 if (overage == 2) // in a different triangle
573 {
574 while (1) {
575 overage = _adjustOverageClassII(&fijk, adjRes, pentLeading4, 1);
576 if (overage != 2) // not in a different triangle
577 break;
578 }
579 }
580
581 // all Class III pentagon edges cross icosa edges
582 // note that Class II pentagons have vertices on the edge,
583 // not edge intersections
584 if (isResClassIII(res) && vert > 0) {
585 // find hex2d of the two vertexes on the last face
586
587 FaceIJK tmpFijk = fijk;
588
589 Vec2d orig2d0;
590 _ijkToHex2d(&lastFijk.coord, &orig2d0);
591
592 int currentToLastDir = adjacentFaceDir[tmpFijk.face][lastFijk.face];
593
594 const FaceOrientIJK* fijkOrient =
595 &faceNeighbors[tmpFijk.face][currentToLastDir];
596
597 tmpFijk.face = fijkOrient->face;
598 CoordIJK* ijk = &tmpFijk.coord;
599
600 // rotate and translate for adjacent face
601 for (int i = 0; i < fijkOrient->ccwRot60; i++) _ijkRotate60ccw(ijk);
602
603 CoordIJK transVec = fijkOrient->translate;
604 _ijkScale(&transVec, unitScaleByCIIres[adjRes] * 3);
605 _ijkAdd(ijk, &transVec, ijk);
606 _ijkNormalize(ijk);
607
608 Vec2d orig2d1;
609 _ijkToHex2d(ijk, &orig2d1);
610
611 // find the appropriate icosa face edge vertexes
612 int maxDim = maxDimByCIIres[adjRes];
613 Vec2d v0 = {3.0 * maxDim, 0.0};
614 Vec2d v1 = {-1.5 * maxDim, 3.0 * M_SQRT3_2 * maxDim};
615 Vec2d v2 = {-1.5 * maxDim, -3.0 * M_SQRT3_2 * maxDim};
616
617 Vec2d* edge0;
618 Vec2d* edge1;
619 switch (adjacentFaceDir[tmpFijk.face][fijk.face]) {
620 case IJ:
621 edge0 = &v0;
622 edge1 = &v1;
623 break;
624 case JK:
625 edge0 = &v1;
626 edge1 = &v2;
627 break;
628 case KI:
629 default:
630 assert(adjacentFaceDir[tmpFijk.face][fijk.face] == KI);
631 edge0 = &v2;
632 edge1 = &v0;
633 break;
634 }
635
636 // find the intersection and add the lat/lon point to the result
637 Vec2d inter;
638 _v2dIntersect(&orig2d0, &orig2d1, edge0, edge1, &inter);
639 _hex2dToGeo(&inter, tmpFijk.face, adjRes, 1,
640 &g->verts[g->numVerts]);
641 g->numVerts++;
642 }
643
644 // convert vertex to lat/lon and add to the result
645 // vert == NUM_PENT_VERTS is only used to test for possible intersection
646 // on last edge
647 if (vert < NUM_PENT_VERTS) {
648 Vec2d vec;
649 _ijkToHex2d(&fijk.coord, &vec);
650 _hex2dToGeo(&vec, fijk.face, adjRes, 1, &g->verts[g->numVerts]);
651 g->numVerts++;
652 }
653
654 lastFijk = fijk;
655 }
656}
657
658/**
659 * Generates the cell boundary in spherical coordinates for a cell given by a
660 * FaceIJK address at a specified resolution.
661 *
662 * @param h The FaceIJK address of the cell.
663 * @param res The H3 resolution of the cell.
664 * @param isPentagon Whether or not the cell is a pentagon.
665 * @param g The spherical coordinates of the cell boundary.
666 */
667void _faceIjkToGeoBoundary(const FaceIJK* h, int res, int isPentagon,
668 GeoBoundary* g) {
669 if (isPentagon) {
670 _faceIjkPentToGeoBoundary(h, res, g);
671 return;
672 }
673
674 // the vertexes of an origin-centered cell in a Class II resolution on a
675 // substrate grid with aperture sequence 33r. The aperture 3 gets us the
676 // vertices, and the 3r gets us back to Class II.
677 // vertices listed ccw from the i-axes
678 CoordIJK vertsCII[NUM_HEX_VERTS] = {
679 {2, 1, 0}, // 0
680 {1, 2, 0}, // 1
681 {0, 2, 1}, // 2
682 {0, 1, 2}, // 3
683 {1, 0, 2}, // 4
684 {2, 0, 1} // 5
685 };
686
687 // the vertexes of an origin-centered cell in a Class III resolution on a
688 // substrate grid with aperture sequence 33r7r. The aperture 3 gets us the
689 // vertices, and the 3r7r gets us to Class II.
690 // vertices listed ccw from the i-axes
691 CoordIJK vertsCIII[NUM_HEX_VERTS] = {
692 {5, 4, 0}, // 0
693 {1, 5, 0}, // 1
694 {0, 5, 4}, // 2
695 {0, 1, 5}, // 3
696 {4, 0, 5}, // 4
697 {5, 0, 1} // 5
698 };
699
700 // get the correct set of substrate vertices for this resolution
701 CoordIJK* verts;
702 if (isResClassIII(res))
703 verts = vertsCIII;
704 else
705 verts = vertsCII;
706
707 // adjust the center point to be in an aperture 33r substrate grid
708 // these should be composed for speed
709 FaceIJK centerIJK = *h;
710 _downAp3(&centerIJK.coord);
711 _downAp3r(&centerIJK.coord);
712
713 // if res is Class III we need to add a cw aperture 7 to get to
714 // icosahedral Class II
715 int adjRes = res;
716 if (isResClassIII(res)) {
717 _downAp7r(&centerIJK.coord);
718 adjRes++;
719 }
720
721 // The center point is now in the same substrate grid as the origin
722 // cell vertices. Add the center point substate coordinates
723 // to each vertex to translate the vertices to that cell.
724 FaceIJK fijkVerts[NUM_HEX_VERTS];
725 for (int v = 0; v < NUM_HEX_VERTS; v++) {
726 fijkVerts[v].face = centerIJK.face;
727 _ijkAdd(&centerIJK.coord, &verts[v], &fijkVerts[v].coord);
728 _ijkNormalize(&fijkVerts[v].coord);
729 }
730
731 // convert each vertex to lat/lon
732 // adjust the face of each vertex as appropriate and introduce
733 // edge-crossing vertices as needed
734 g->numVerts = 0;
735 int lastFace = -1;
736 int lastOverage = 0; // 0: none; 1: edge; 2: overage
737 for (int vert = 0; vert < NUM_HEX_VERTS + 1; vert++) {
738 int v = vert % NUM_HEX_VERTS;
739
740 FaceIJK fijk = fijkVerts[v];
741
742 int pentLeading4 = 0;
743 int overage = _adjustOverageClassII(&fijk, adjRes, pentLeading4, 1);
744
745 /*
746 Check for edge-crossing. Each face of the underlying icosahedron is a
747 different projection plane. So if an edge of the hexagon crosses an
748 icosahedron edge, an additional vertex must be introduced at that
749 intersection point. Then each half of the cell edge can be projected
750 to geographic coordinates using the appropriate icosahedron face
751 projection. Note that Class II cell edges have vertices on the face
752 edge, with no edge line intersections.
753 */
754 if (isResClassIII(res) && vert > 0 && fijk.face != lastFace &&
755 lastOverage != 1) {
756 // find hex2d of the two vertexes on original face
757 int lastV = (v + 5) % NUM_HEX_VERTS;
758 Vec2d orig2d0;
759 _ijkToHex2d(&fijkVerts[lastV].coord, &orig2d0);
760
761 Vec2d orig2d1;
762 _ijkToHex2d(&fijkVerts[v].coord, &orig2d1);
763
764 // find the appropriate icosa face edge vertexes
765 int maxDim = maxDimByCIIres[adjRes];
766 Vec2d v0 = {3.0 * maxDim, 0.0};
767 Vec2d v1 = {-1.5 * maxDim, 3.0 * M_SQRT3_2 * maxDim};
768 Vec2d v2 = {-1.5 * maxDim, -3.0 * M_SQRT3_2 * maxDim};
769
770 int face2 = ((lastFace == centerIJK.face) ? fijk.face : lastFace);
771 Vec2d* edge0;
772 Vec2d* edge1;
773 switch (adjacentFaceDir[centerIJK.face][face2]) {
774 case IJ:
775 edge0 = &v0;
776 edge1 = &v1;
777 break;
778 case JK:
779 edge0 = &v1;
780 edge1 = &v2;
781 break;
782 case KI:
783 default:
784 assert(adjacentFaceDir[centerIJK.face][face2] == KI);
785 edge0 = &v2;
786 edge1 = &v0;
787 break;
788 }
789
790 // find the intersection and add the lat/lon point to the result
791 Vec2d inter;
792 _v2dIntersect(&orig2d0, &orig2d1, edge0, edge1, &inter);
793 /*
794 If a point of intersection occurs at a hexagon vertex, then each
795 adjacent hexagon edge will lie completely on a single icosahedron
796 face, and no additional vertex is required.
797 */
798 bool isIntersectionAtVertex =
799 _v2dEquals(&orig2d0, &inter) || _v2dEquals(&orig2d1, &inter);
800 if (!isIntersectionAtVertex) {
801 _hex2dToGeo(&inter, centerIJK.face, adjRes, 1,
802 &g->verts[g->numVerts]);
803 g->numVerts++;
804 }
805 }
806
807 // convert vertex to lat/lon and add to the result
808 // vert == NUM_HEX_VERTS is only used to test for possible intersection
809 // on last edge
810 if (vert < NUM_HEX_VERTS) {
811 Vec2d vec;
812 _ijkToHex2d(&fijk.coord, &vec);
813 _hex2dToGeo(&vec, fijk.face, adjRes, 1, &g->verts[g->numVerts]);
814 g->numVerts++;
815 }
816
817 lastFace = fijk.face;
818 lastOverage = overage;
819 }
820}
821
822/**
823 * Adjusts a FaceIJK address in place so that the resulting cell address is
824 * relative to the correct icosahedral face.
825 *
826 * @param fijk The FaceIJK address of the cell.
827 * @param res The H3 resolution of the cell.
828 * @param pentLeading4 Whether or not the cell is a pentagon with a leading
829 * digit 4.
830 * @param substrate Whether or not the cell is in a substrate grid.
831 * @return 0 if on original face (no overage); 1 if on face edge (only occurs
832 * on substrate grids); 2 if overage on new face interior
833 */
834int _adjustOverageClassII(FaceIJK* fijk, int res, int pentLeading4,
835 int substrate) {
836 int overage = 0;
837
838 CoordIJK* ijk = &fijk->coord;
839
840 // get the maximum dimension value; scale if a substrate grid
841 int maxDim = maxDimByCIIres[res];
842 if (substrate) maxDim *= 3;
843
844 // check for overage
845 if (substrate && ijk->i + ijk->j + ijk->k == maxDim) // on edge
846 overage = 1;
847 else if (ijk->i + ijk->j + ijk->k > maxDim) // overage
848 {
849 overage = 2;
850
851 const FaceOrientIJK* fijkOrient;
852 if (ijk->k > 0) {
853 if (ijk->j > 0) // jk "quadrant"
854 fijkOrient = &faceNeighbors[fijk->face][JK];
855 else // ik "quadrant"
856 {
857 fijkOrient = &faceNeighbors[fijk->face][KI];
858
859 // adjust for the pentagonal missing sequence
860 if (pentLeading4) {
861 // translate origin to center of pentagon
862 CoordIJK origin;
863 _setIJK(&origin, maxDim, 0, 0);
864 CoordIJK tmp;
865 _ijkSub(ijk, &origin, &tmp);
866 // rotate to adjust for the missing sequence
867 _ijkRotate60cw(&tmp);
868 // translate the origin back to the center of the triangle
869 _ijkAdd(&tmp, &origin, ijk);
870 }
871 }
872 } else // ij "quadrant"
873 fijkOrient = &faceNeighbors[fijk->face][IJ];
874
875 fijk->face = fijkOrient->face;
876
877 // rotate and translate for adjacent face
878 for (int i = 0; i < fijkOrient->ccwRot60; i++) _ijkRotate60ccw(ijk);
879
880 CoordIJK transVec = fijkOrient->translate;
881 int unitScale = unitScaleByCIIres[res];
882 if (substrate) unitScale *= 3;
883 _ijkScale(&transVec, unitScale);
884 _ijkAdd(ijk, &transVec, ijk);
885 _ijkNormalize(ijk);
886
887 // overage points on pentagon boundaries can end up on edges
888 if (substrate && ijk->i + ijk->j + ijk->k == maxDim) // on edge
889 overage = 1;
890 }
891
892 return overage;
893}
894