/*

 * Copyright 2016-2023, 2026 Uber Technologies, Inc.

 *

 * Licensed under the Apache License, Version 2.0 (the "License");

 * you may not use this file except in compliance with the License.

 * You may obtain a copy of the License at

 *

 *         http://www.apache.org/licenses/LICENSE-2.0

 *

 * Unless required by applicable law or agreed to in writing, software

 * distributed under the License is distributed on an "AS IS" BASIS,

 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.

 * See the License for the specific language governing permissions and

 * limitations under the License.

 */

/** @file faceijk.c

 * @brief   Functions for working with icosahedral face-centered hex IJK

 *  coordinate systems.

 */

#include "faceijk.h"

#include <assert.h>

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "constants.h"

#include "coordijk.h"

#include "h3Index.h"

#include "latLng.h"

#include "vec3d.h"

/** square root of 7 and inverse square root of 7 */

#define M_SQRT7 2.6457513110645905905016157536392604257102

#define M_RSQRT7 0.37796447300922722721451653623418006081576

/** @brief icosahedron face centers in x/y/z on the unit sphere */

static const Vec3d faceCenterPoint[NUM_ICOSA_FACES] = {

    {0.2199307791404606, 0.6583691780274996, 0.7198475378926182},     // face  0

    {-0.2139234834501421, 0.1478171829550703, 0.9656017935214205},    // face  1

    {0.1092625278784797, -0.4811951572873210, 0.8697775121287253},    // face  2

    {0.7428567301586791, -0.3593941678278028, 0.5648005936517033},    // face  3

    {0.8112534709140969, 0.3448953237639384, 0.4721387736413930},     // face  4

    {-0.1055498149613921, 0.9794457296411413, 0.1718874610009365},    // face  5

    {-0.8075407579970092, 0.1533552485898818, 0.5695261994882688},    // face  6

    {-0.2846148069787907, -0.8644080972654206, 0.4144792552473539},   // face  7

    {0.7405621473854482, -0.6673299564565524, -0.0789837646326737},   // face  8

    {0.8512303986474293, 0.4722343788582681, -0.2289137388687808},    // face  9

    {-0.7405621473854481, 0.6673299564565524, 0.0789837646326737},    // face 10

    {-0.8512303986474292, -0.4722343788582682, 0.2289137388687808},   // face 11

    {0.1055498149613919, -0.9794457296411413, -0.1718874610009365},   // face 12

    {0.8075407579970092, -0.1533552485898819, -0.5695261994882688},   // face 13

    {0.2846148069787908, 0.8644080972654204, -0.4144792552473539},    // face 14

    {-0.7428567301586791, 0.3593941678278027, -0.5648005936517033},   // face 15

    {-0.8112534709140971, -0.3448953237639382, -0.4721387736413930},  // face 16

    {-0.2199307791404607, -0.6583691780274996, -0.7198475378926182},  // face 17

    {0.2139234834501420, -0.1478171829550704, -0.9656017935214205},   // face 18

    {-0.1092625278784796, 0.4811951572873210, -0.8697775121287253},   // face 19

};

/** @brief icosahedron face ijk axes as azimuth in radians from face center to

 * vertex 0/1/2 respectively

 */

static const double faceAxesAzRadsCII[NUM_ICOSA_FACES][3] = {

    {5.619958268523939882, 3.525563166130744542,

     1.431168063737548730},  // face  0

    {5.760339081714187279, 3.665943979320991689,

     1.571548876927796127},  // face  1

    {0.780213654393430055, 4.969003859179821079,

     2.874608756786625655},  // face  2

    {0.430469363979999913, 4.619259568766391033,

     2.524864466373195467},  // face  3

    {6.130269123335111400, 4.035874020941915804,

     1.941478918548720291},  // face  4

    {2.692877706530642877, 0.598482604137447119,

     4.787272808923838195},  // face  5

    {2.982963003477243874, 0.888567901084048369,

     5.077358105870439581},  // face  6

    {3.532912002790141181, 1.438516900396945656,

     5.627307105183336758},  // face  7

    {3.494305004259568154, 1.399909901866372864,

     5.588700106652763840},  // face  8

    {3.003214169499538391, 0.908819067106342928,

     5.097609271892733906},  // face  9

    {5.930472956509811562, 3.836077854116615875,

     1.741682751723420374},  // face 10

    {0.138378484090254847, 4.327168688876645809,

     2.232773586483450311},  // face 11

    {0.448714947059150361, 4.637505151845541521,

     2.543110049452346120},  // face 12

    {0.158629650112549365, 4.347419854898940135,

     2.253024752505744869},  // face 13

    {5.891865957979238535, 3.797470855586042958,

     1.703075753192847583},  // face 14

    {2.711123289609793325, 0.616728187216597771,

     4.805518392002988683},  // face 15

    {3.294508837434268316, 1.200113735041072948,

     5.388903939827463911},  // face 16

    {3.804819692245439833, 1.710424589852244509,

     5.899214794638635174},  // face 17

    {3.664438879055192436, 1.570043776661997111,

     5.758833981448388027},  // face 18

    {2.361378999196363184, 0.266983896803167583,

     4.455774101589558636},  // face 19

};

/** @brief Definition of which faces neighbor each other. */

static const FaceOrientIJK faceNeighbors[NUM_ICOSA_FACES][4] = {

    {

        // face 0

        {0, {0, 0, 0}, 0},  // central face

        {4, {2, 0, 2}, 1},  // ij quadrant

        {1, {2, 2, 0}, 5},  // ki quadrant

        {5, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 1

        {1, {0, 0, 0}, 0},  // central face

        {0, {2, 0, 2}, 1},  // ij quadrant

        {2, {2, 2, 0}, 5},  // ki quadrant

        {6, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 2

        {2, {0, 0, 0}, 0},  // central face

        {1, {2, 0, 2}, 1},  // ij quadrant

        {3, {2, 2, 0}, 5},  // ki quadrant

        {7, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 3

        {3, {0, 0, 0}, 0},  // central face

        {2, {2, 0, 2}, 1},  // ij quadrant

        {4, {2, 2, 0}, 5},  // ki quadrant

        {8, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 4

        {4, {0, 0, 0}, 0},  // central face

        {3, {2, 0, 2}, 1},  // ij quadrant

        {0, {2, 2, 0}, 5},  // ki quadrant

        {9, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 5

        {5, {0, 0, 0}, 0},   // central face

        {10, {2, 2, 0}, 3},  // ij quadrant

        {14, {2, 0, 2}, 3},  // ki quadrant

        {0, {0, 2, 2}, 3}    // jk quadrant

    },

    {

        // face 6

        {6, {0, 0, 0}, 0},   // central face

        {11, {2, 2, 0}, 3},  // ij quadrant

        {10, {2, 0, 2}, 3},  // ki quadrant

        {1, {0, 2, 2}, 3}    // jk quadrant

    },

    {

        // face 7

        {7, {0, 0, 0}, 0},   // central face

        {12, {2, 2, 0}, 3},  // ij quadrant

        {11, {2, 0, 2}, 3},  // ki quadrant

        {2, {0, 2, 2}, 3}    // jk quadrant

    },

    {

        // face 8

        {8, {0, 0, 0}, 0},   // central face

        {13, {2, 2, 0}, 3},  // ij quadrant

        {12, {2, 0, 2}, 3},  // ki quadrant

        {3, {0, 2, 2}, 3}    // jk quadrant

    },

    {

        // face 9

        {9, {0, 0, 0}, 0},   // central face

        {14, {2, 2, 0}, 3},  // ij quadrant

        {13, {2, 0, 2}, 3},  // ki quadrant

        {4, {0, 2, 2}, 3}    // jk quadrant

    },

    {

        // face 10

        {10, {0, 0, 0}, 0},  // central face

        {5, {2, 2, 0}, 3},   // ij quadrant

        {6, {2, 0, 2}, 3},   // ki quadrant

        {15, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 11

        {11, {0, 0, 0}, 0},  // central face

        {6, {2, 2, 0}, 3},   // ij quadrant

        {7, {2, 0, 2}, 3},   // ki quadrant

        {16, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 12

        {12, {0, 0, 0}, 0},  // central face

        {7, {2, 2, 0}, 3},   // ij quadrant

        {8, {2, 0, 2}, 3},   // ki quadrant

        {17, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 13

        {13, {0, 0, 0}, 0},  // central face

        {8, {2, 2, 0}, 3},   // ij quadrant

        {9, {2, 0, 2}, 3},   // ki quadrant

        {18, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 14

        {14, {0, 0, 0}, 0},  // central face

        {9, {2, 2, 0}, 3},   // ij quadrant

        {5, {2, 0, 2}, 3},   // ki quadrant

        {19, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 15

        {15, {0, 0, 0}, 0},  // central face

        {16, {2, 0, 2}, 1},  // ij quadrant

        {19, {2, 2, 0}, 5},  // ki quadrant

        {10, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 16

        {16, {0, 0, 0}, 0},  // central face

        {17, {2, 0, 2}, 1},  // ij quadrant

        {15, {2, 2, 0}, 5},  // ki quadrant

        {11, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 17

        {17, {0, 0, 0}, 0},  // central face

        {18, {2, 0, 2}, 1},  // ij quadrant

        {16, {2, 2, 0}, 5},  // ki quadrant

        {12, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 18

        {18, {0, 0, 0}, 0},  // central face

        {19, {2, 0, 2}, 1},  // ij quadrant

        {17, {2, 2, 0}, 5},  // ki quadrant

        {13, {0, 2, 2}, 3}   // jk quadrant

    },

    {

        // face 19

        {19, {0, 0, 0}, 0},  // central face

        {15, {2, 0, 2}, 1},  // ij quadrant

        {18, {2, 2, 0}, 5},  // ki quadrant

        {14, {0, 2, 2}, 3}   // jk quadrant

    }};

/** @brief direction from the origin face to the destination face, relative to

 * the origin face's coordinate system, or -1 if not adjacent.

 */

static const int adjacentFaceDir[NUM_ICOSA_FACES][NUM_ICOSA_FACES] = {

    {0,  KI, -1, -1, IJ, JK, -1, -1, -1, -1,

     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},  // face 0

    {IJ, 0,  KI, -1, -1, -1, JK, -1, -1, -1,

     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},  // face 1

    {-1, IJ, 0,  KI, -1, -1, -1, JK, -1, -1,

     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},  // face 2

    {-1, -1, IJ, 0,  KI, -1, -1, -1, JK, -1,

     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},  // face 3

    {KI, -1, -1, IJ, 0,  -1, -1, -1, -1, JK,

     -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},  // face 4

    {JK, -1, -1, -1, -1, 0,  -1, -1, -1, -1,

     IJ, -1, -1, -1, KI, -1, -1, -1, -1, -1},  // face 5

    {-1, JK, -1, -1, -1, -1, 0,  -1, -1, -1,

     KI, IJ, -1, -1, -1, -1, -1, -1, -1, -1},  // face 6

    {-1, -1, JK, -1, -1, -1, -1, 0,  -1, -1,

     -1, KI, IJ, -1, -1, -1, -1, -1, -1, -1},  // face 7

    {-1, -1, -1, JK, -1, -1, -1, -1, 0,  -1,

     -1, -1, KI, IJ, -1, -1, -1, -1, -1, -1},  // face 8

    {-1, -1, -1, -1, JK, -1, -1, -1, -1, 0,

     -1, -1, -1, KI, IJ, -1, -1, -1, -1, -1},  // face 9

    {-1, -1, -1, -1, -1, IJ, KI, -1, -1, -1,

     0,  -1, -1, -1, -1, JK, -1, -1, -1, -1},  // face 10

    {-1, -1, -1, -1, -1, -1, IJ, KI, -1, -1,

     -1, 0,  -1, -1, -1, -1, JK, -1, -1, -1},  // face 11

    {-1, -1, -1, -1, -1, -1, -1, IJ, KI, -1,

     -1, -1, 0,  -1, -1, -1, -1, JK, -1, -1},  // face 12

    {-1, -1, -1, -1, -1, -1, -1, -1, IJ, KI,

     -1, -1, -1, 0,  -1, -1, -1, -1, JK, -1},  // face 13

    {-1, -1, -1, -1, -1, KI, -1, -1, -1, IJ,

     -1, -1, -1, -1, 0,  -1, -1, -1, -1, JK},  // face 14

    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,

     JK, -1, -1, -1, -1, 0,  IJ, -1, -1, KI},  // face 15

    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,

     -1, JK, -1, -1, -1, KI, 0,  IJ, -1, -1},  // face 16

    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,

     -1, -1, JK, -1, -1, -1, KI, 0,  IJ, -1},  // face 17

    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,

     -1, -1, -1, JK, -1, -1, -1, KI, 0,  IJ},  // face 18

    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,

     -1, -1, -1, -1, JK, IJ, -1, -1, KI, 0}  // face 19

};

/** @brief overage distance table */

static const int maxDimByCIIres[] = {

    2,        // res  0

    -1,       // res  1

    14,       // res  2

    -1,       // res  3

    98,       // res  4

    -1,       // res  5

    686,      // res  6

    -1,       // res  7

    4802,     // res  8

    -1,       // res  9

    33614,    // res 10

    -1,       // res 11

    235298,   // res 12

    -1,       // res 13

    1647086,  // res 14

    -1,       // res 15

    11529602  // res 16

};

/** @brief unit scale distance table */

static const int unitScaleByCIIres[] = {

    1,       // res  0

    -1,      // res  1

    7,       // res  2

    -1,      // res  3

    49,      // res  4

    -1,      // res  5

    343,     // res  6

    -1,      // res  7

    2401,    // res  8

    -1,      // res  9

    16807,   // res 10

    -1,      // res 11

    117649,  // res 12

    -1,      // res 13

    823543,  // res 14

    -1,      // res 15

    5764801  // res 16

};

// Forward declares to make diff nicer

// TODO: remove and reorder functions after landing

static void _vec3ToHex2d(const Vec3d *p, int res, int *face, Vec2d *v);

static void _vec3ToClosestFace(const Vec3d *v, int *face, double *sqd);

/**

 * Encodes a Vec3d coordinate to the FaceIJK address of the containing

 * cell at the specified resolution.

 *

 * Vec3d p is expected to be on the unit sphere.

 *

 * @param p The Vec3d coordinates to encode.

 * @param res The desired H3 resolution for the encoding.

 * @param h Output: FaceIJK address of the containing cell at resolution res.

 */

void _vec3ToFaceIjk(Vec3d p, int res, FaceIJK *h) {

    // first convert to hex2d

    Vec2d v;

    _vec3ToHex2d(&p, res, &h->face, &v);

    // then convert to ijk+

    _hex2dToCoordIJK(&v, &h->coord);

}

/**

 * Compute the local north and east directions on the tangent plane

 * at a point on the unit sphere.

 *

 * Will not work if p is at a pole, but icosahedron face centers

 * are never at the poles.

 *

 * @param p Unit vector on the sphere.

 * @param north Output: local north direction on tangent plane.

 * @param east Output: local east direction on tangent plane.

 */

static inline void _vec3TangentBasis(Vec3d p, Vec3d *north, Vec3d *east) {

    Vec3d northPole = {0.0, 0.0, 1.0};

    *north = vec3LinComb(1.0, northPole, -vec3Dot(northPole, p), p);

    vec3Normalize(north);

    *east = vec3Cross(*north, p);

}

/**

 * Calculates the azimuth from p1 to p2.

 * @param p1 The first vector.

 * @param p2 The second vector.

 * @return The azimuth in radians.

 */

static inline double _vec3AzimuthRads(Vec3d p1, Vec3d p2) {

    Vec3d northDir, eastDir;

    _vec3TangentBasis(p1, &northDir, &eastDir);

    // project p2 onto tangent plane at p1

    Vec3d p2Proj = vec3LinComb(1.0, p2, -vec3Dot(p2, p1), p1);

    vec3Normalize(&p2Proj);

    return atan2(vec3Dot(p2Proj, eastDir), vec3Dot(p2Proj, northDir));

}

/**

 * Encodes a coordinate on the sphere to the corresponding icosahedral face and

 * containing 2D hex coordinates relative to that face center.

 *

 * Vec3d p is expected to be on the unit sphere.

 *

 * @param p The Vec3d coordinates to encode.

 * @param res The desired H3 resolution for the encoding.

 * @param face Output: The icosahedral face containing the coordinates.

 * @param v Output: The 2D hex coordinates of the cell containing the point.

 */

static void _vec3ToHex2d(const Vec3d *p, int res, int *face, Vec2d *v) {

    // determine the icosahedron face

    double sqd;

    _vec3ToClosestFace(p, face, &sqd);

    // cos(r) = 1 - 2 * sin^2(r/2) = 1 - 2 * (sqd / 4) = 1 - sqd/2

    double r = acos(1 - sqd * 0.5);

    if (r < EPSILON) {

        v->x = v->y = 0.0;

        return;

    }

    // now have face and r, now find CCW theta from CII i-axis

    double theta = _posAngleRads(

        faceAxesAzRadsCII[*face][0] -

        _posAngleRads(_vec3AzimuthRads(faceCenterPoint[*face], *p)));

    // adjust theta for Class III (odd resolutions)

    if (isResolutionClassIII(res))

        theta = _posAngleRads(theta - M_AP7_ROT_RADS);

    // perform gnomonic scaling of r

    r = tan(r);

    // scale for current resolution length u

    r *= INV_RES0_U_GNOMONIC;

    for (int i = 0; i < res; i++) r *= M_SQRT7;

    // we now have (r, theta) in hex2d with theta ccw from x-axes

    // convert to local x,y

    v->x = r * cos(theta);

    v->y = r * sin(theta);

}

/**

 * Determines the 3D coordinates of a cell given by 2D

 * hex coordinates on a particular icosahedral face.

 *

 * @param v The 2D hex coordinates of the cell.

 * @param face The icosahedral face upon which the 2D hex coordinate system is

 *             centered.

 * @param res The H3 resolution of the cell.

 * @param substrate Indicates whether or not this grid is actually a substrate

 *        grid relative to the specified resolution.

 * @param v3 Output: the 3D coordinates of the cell center point

 */

static void _hex2dToVec3(const Vec2d *v, int face, int res, int substrate,

                         Vec3d *v3) {

    // calculate (r, theta) in hex2d

    double r = _v2dMag(v);

    if (r < EPSILON) {

        *v3 = faceCenterPoint[face];

        return;

    }

    double theta = atan2(v->y, v->x);

    // scale for current resolution length u

    for (int i = 0; i < res; i++) r *= M_RSQRT7;

    // scale accordingly if this is a substrate grid

    if (substrate) {

        r *= M_ONETHIRD;

        if (isResolutionClassIII(res)) r *= M_RSQRT7;

    }

    r *= RES0_U_GNOMONIC;

    // perform inverse gnomonic scaling of r

    r = atan(r);

    // adjust theta for Class III

    // if a substrate grid, then it's already been adjusted for Class III

    if (!substrate && isResolutionClassIII(res))

        theta = _posAngleRads(theta + M_AP7_ROT_RADS);

    // find theta as an azimuth

    theta = _posAngleRads(faceAxesAzRadsCII[face][0] - theta);

    // now find the point at (r,theta) from the face center

    Vec3d northDir, eastDir;

    _vec3TangentBasis(faceCenterPoint[face], &northDir, &eastDir);

    Vec3d dir = vec3LinComb(cos(theta), northDir, sin(theta), eastDir);

    *v3 = vec3LinComb(cos(r), faceCenterPoint[face], sin(r), dir);

    vec3Normalize(v3);

}

/**

 * Determines the center point in 3D coordinates of a cell given by

 * a FaceIJK address at a specified resolution.

 *

 * @param h The FaceIJK address of the cell.

 * @param res The H3 resolution of the cell.

 * @param g Output: The 3D coordinates of the cell center point.

 */

void _faceIjkToVec3(const FaceIJK *h, int res, Vec3d *g) {

    Vec2d v;

    _ijkToHex2d(&h->coord, &v);

    _hex2dToVec3(&v, h->face, res, 0, g);

}

/**

 * Generates the cell boundary in spherical coordinates for a pentagonal cell

 * given by a FaceIJK address at a specified resolution.

 *

 * @param h The FaceIJK address of the pentagonal cell.

 * @param res The H3 resolution of the cell.

 * @param start The first topological vertex to return.

 * @param length The number of topological vertexes to return.

 * @param g Output: The spherical coordinates of the cell boundary.

 */

void _faceIjkPentToCellBoundary(const FaceIJK *h, int res, int start,

                                int length, CellBoundary *g) {

    int adjRes = res;

    FaceIJK centerIJK = *h;

    FaceIJK fijkVerts[NUM_PENT_VERTS];

    _faceIjkPentToVerts(&centerIJK, &adjRes, fijkVerts);

    // If we're returning the entire loop, we need one more iteration in case

    // of a distortion vertex on the last edge

    int additionalIteration = length == NUM_PENT_VERTS ? 1 : 0;

    // convert each vertex to lat/lng

    // adjust the face of each vertex as appropriate and introduce

    // edge-crossing vertices as needed

    g->numVerts = 0;

    FaceIJK lastFijk = {0};

    for (int vert = start; vert < start + length + additionalIteration;

         vert++) {

        int v = vert % NUM_PENT_VERTS;

        FaceIJK fijk = fijkVerts[v];

        _adjustPentVertOverage(&fijk, adjRes);

        // all Class III pentagon edges cross icosa edges

        // note that Class II pentagons have vertices on the edge,

        // not edge intersections

        if (isResolutionClassIII(res) && vert > start) {

            // find hex2d of the two vertexes on the last face

            FaceIJK tmpFijk = fijk;

            Vec2d orig2d0;

            _ijkToHex2d(&lastFijk.coord, &orig2d0);

            int currentToLastDir = adjacentFaceDir[tmpFijk.face][lastFijk.face];

            const FaceOrientIJK *fijkOrient =

                &faceNeighbors[tmpFijk.face][currentToLastDir];

            tmpFijk.face = fijkOrient->face;

            CoordIJK *ijk = &tmpFijk.coord;

            // rotate and translate for adjacent face

            for (int i = 0; i < fijkOrient->ccwRot60; i++) _ijkRotate60ccw(ijk);

            CoordIJK transVec = fijkOrient->translate;

            _ijkScale(&transVec, unitScaleByCIIres[adjRes] * 3);

            _ijkAdd(ijk, &transVec, ijk);

            _ijkNormalize(ijk);

            Vec2d orig2d1;

            _ijkToHex2d(ijk, &orig2d1);

            // find the appropriate icosa face edge vertexes

            int maxDim = maxDimByCIIres[adjRes];

            Vec2d v0 = {3.0 * maxDim, 0.0};

            Vec2d v1 = {-1.5 * maxDim, 3.0 * M_SQRT3_2 * maxDim};

            Vec2d v2 = {-1.5 * maxDim, -3.0 * M_SQRT3_2 * maxDim};

            Vec2d *edge0;

            Vec2d *edge1;

            switch (adjacentFaceDir[tmpFijk.face][fijk.face]) {

                case IJ:

                    edge0 = &v0;

                    edge1 = &v1;

                    break;

                case JK:

                    edge0 = &v1;

                    edge1 = &v2;

                    break;

                case KI:

                default:

                    assert(adjacentFaceDir[tmpFijk.face][fijk.face] == KI);

                    edge0 = &v2;

                    edge1 = &v0;

                    break;

            }

            // find the intersection and add the lat/lng point to the result

            Vec2d inter;

            _v2dIntersect(&orig2d0, &orig2d1, edge0, edge1, &inter);

            Vec3d v3;

            _hex2dToVec3(&inter, tmpFijk.face, adjRes, 1, &v3);

            g->verts[g->numVerts] = vec3ToLatLng(v3);

            g->numVerts++;

        }

        // convert vertex to lat/lng and add to the result

        // vert == start + NUM_PENT_VERTS is only used to test for possible

        // intersection on last edge

        if (vert < start + NUM_PENT_VERTS) {

            Vec2d vec;

            _ijkToHex2d(&fijk.coord, &vec);

            Vec3d v3;

            _hex2dToVec3(&vec, fijk.face, adjRes, 1, &v3);

            g->verts[g->numVerts] = vec3ToLatLng(v3);

            g->numVerts++;

        }

        lastFijk = fijk;

    }

}

/**

 * Get the vertices of a pentagon cell as substrate FaceIJK addresses

 *

 * @param fijk The FaceIJK address of the cell.

 * @param res In/out: the H3 resolution of the cell, adjusted for substrate.

 * @param fijkVerts Output: array for the vertices.

 */

void _faceIjkPentToVerts(FaceIJK *fijk, int *res, FaceIJK *fijkVerts) {

    // the vertexes of an origin-centered pentagon in a Class II resolution on a

    // substrate grid with aperture sequence 33r. The aperture 3 gets us the

    // vertices, and the 3r gets us back to Class II.

    // vertices listed ccw from the i-axes

    CoordIJK vertsCII[NUM_PENT_VERTS] = {

        {2, 1, 0},  // 0

        {1, 2, 0},  // 1

        {0, 2, 1},  // 2

        {0, 1, 2},  // 3

        {1, 0, 2},  // 4

    };

    // the vertexes of an origin-centered pentagon in a Class III resolution on

    // a substrate grid with aperture sequence 33r7r. The aperture 3 gets us the

    // vertices, and the 3r7r gets us to Class II. vertices listed ccw from the

    // i-axes

    CoordIJK vertsCIII[NUM_PENT_VERTS] = {

        {5, 4, 0},  // 0

        {1, 5, 0},  // 1

        {0, 5, 4},  // 2

        {0, 1, 5},  // 3

        {4, 0, 5},  // 4

    };

    // get the correct set of substrate vertices for this resolution

    CoordIJK *verts;

    if (isResolutionClassIII(*res))

        verts = vertsCIII;

    else

        verts = vertsCII;

    // adjust the center point to be in an aperture 33r substrate grid

    // these should be composed for speed

    _downAp3(&fijk->coord);

    _downAp3r(&fijk->coord);

    // if res is Class III we need to add a cw aperture 7 to get to

    // icosahedral Class II

    if (isResolutionClassIII(*res)) {

        _downAp7r(&fijk->coord);

        *res += 1;

    }

    // The center point is now in the same substrate grid as the origin

    // cell vertices. Add the center point substate coordinates

    // to each vertex to translate the vertices to that cell.

    for (int v = 0; v < NUM_PENT_VERTS; v++) {

        fijkVerts[v].face = fijk->face;

        _ijkAdd(&fijk->coord, &verts[v], &fijkVerts[v].coord);

        _ijkNormalize(&fijkVerts[v].coord);

    }

}

/**

 * Generates the cell boundary in spherical coordinates for a cell given by a

 * FaceIJK address at a specified resolution.

 *

 * @param h The FaceIJK address of the cell.

 * @param res The H3 resolution of the cell.

 * @param start The first topological vertex to return.

 * @param length The number of topological vertexes to return.

 * @param g Output: The spherical coordinates of the cell boundary.

 */

void _faceIjkToCellBoundary(const FaceIJK *h, int res, int start, int length,

                            CellBoundary *g) {

    int adjRes = res;

    FaceIJK centerIJK = *h;

    FaceIJK fijkVerts[NUM_HEX_VERTS];

    _faceIjkToVerts(&centerIJK, &adjRes, fijkVerts);

    // If we're returning the entire loop, we need one more iteration in case

    // of a distortion vertex on the last edge

    int additionalIteration = length == NUM_HEX_VERTS ? 1 : 0;

    // convert each vertex to lat/lng

    // adjust the face of each vertex as appropriate and introduce

    // edge-crossing vertices as needed

    g->numVerts = 0;

    int lastFace = -1;

    Overage lastOverage = NO_OVERAGE;

    for (int vert = start; vert < start + length + additionalIteration;

         vert++) {

        int v = vert % NUM_HEX_VERTS;

        FaceIJK fijk = fijkVerts[v];

        const int pentLeading4 = 0;

        Overage overage = _adjustOverageClassII(&fijk, adjRes, pentLeading4, 1);

        /*

        Check for edge-crossing. Each face of the underlying icosahedron is a

        different projection plane. So if an edge of the hexagon crosses an

        icosahedron edge, an additional vertex must be introduced at that

        intersection point. Then each half of the cell edge can be projected

        to geographic coordinates using the appropriate icosahedron face

        projection. Note that Class II cell edges have vertices on the face

        edge, with no edge line intersections.

        */

        if (isResolutionClassIII(res) && vert > start &&

            fijk.face != lastFace && lastOverage != FACE_EDGE) {

            // find hex2d of the two vertexes on original face

            int lastV = (v + 5) % NUM_HEX_VERTS;

            Vec2d orig2d0;

            _ijkToHex2d(&fijkVerts[lastV].coord, &orig2d0);

            Vec2d orig2d1;

            _ijkToHex2d(&fijkVerts[v].coord, &orig2d1);

            // find the appropriate icosa face edge vertexes

            int maxDim = maxDimByCIIres[adjRes];

            Vec2d v0 = {3.0 * maxDim, 0.0};

            Vec2d v1 = {-1.5 * maxDim, 3.0 * M_SQRT3_2 * maxDim};

            Vec2d v2 = {-1.5 * maxDim, -3.0 * M_SQRT3_2 * maxDim};

            int face2 = ((lastFace == centerIJK.face) ? fijk.face : lastFace);

            Vec2d *edge0;

            Vec2d *edge1;

            switch (adjacentFaceDir[centerIJK.face][face2]) {

                case IJ:

                    edge0 = &v0;

                    edge1 = &v1;

                    break;

                case JK:

                    edge0 = &v1;

                    edge1 = &v2;

                    break;

                // case KI:

                default:

                    assert(adjacentFaceDir[centerIJK.face][face2] == KI);

                    edge0 = &v2;

                    edge1 = &v0;

                    break;

            }

            // find the intersection and add the lat/lng point to the result

            Vec2d inter;

            _v2dIntersect(&orig2d0, &orig2d1, edge0, edge1, &inter);

            /*

            If a point of intersection occurs at a hexagon vertex, then each

            adjacent hexagon edge will lie completely on a single icosahedron

            face, and no additional vertex is required.

            */

            bool isIntersectionAtVertex = _v2dAlmostEquals(&orig2d0, &inter) ||

                                          _v2dAlmostEquals(&orig2d1, &inter);

            if (!isIntersectionAtVertex) {

                Vec3d v3;

                _hex2dToVec3(&inter, centerIJK.face, adjRes, 1, &v3);

                g->verts[g->numVerts] = vec3ToLatLng(v3);

                g->numVerts++;

            }

        }

        // convert vertex to lat/lng and add to the result

        // vert == start + NUM_HEX_VERTS is only used to test for possible

        // intersection on last edge

        if (vert < start + NUM_HEX_VERTS) {

            Vec2d vec;

            _ijkToHex2d(&fijk.coord, &vec);

            Vec3d v3;

            _hex2dToVec3(&vec, fijk.face, adjRes, 1, &v3);

            g->verts[g->numVerts] = vec3ToLatLng(v3);

            g->numVerts++;

        }

        lastFace = fijk.face;

        lastOverage = overage;

    }

}

/**

 * Get the vertices of a cell as substrate FaceIJK addresses

 *

 * @param fijk The FaceIJK address of the cell.

 * @param res In/out: the H3 resolution of the cell, adjusted for substrate.

 * @param fijkVerts Output: array for the vertices.

 */

void _faceIjkToVerts(FaceIJK *fijk, int *res, FaceIJK *fijkVerts) {

    // the vertexes of an origin-centered cell in a Class II resolution on a

    // substrate grid with aperture sequence 33r. The aperture 3 gets us the

    // vertices, and the 3r gets us back to Class II.

    // vertices listed ccw from the i-axes

    CoordIJK vertsCII[NUM_HEX_VERTS] = {

        {2, 1, 0},  // 0

        {1, 2, 0},  // 1

        {0, 2, 1},  // 2

        {0, 1, 2},  // 3

        {1, 0, 2},  // 4

        {2, 0, 1}   // 5

    };

    // the vertexes of an origin-centered cell in a Class III resolution on a

    // substrate grid with aperture sequence 33r7r. The aperture 3 gets us the

    // vertices, and the 3r7r gets us to Class II.

    // vertices listed ccw from the i-axes

    CoordIJK vertsCIII[NUM_HEX_VERTS] = {

        {5, 4, 0},  // 0

        {1, 5, 0},  // 1

        {0, 5, 4},  // 2

        {0, 1, 5},  // 3

        {4, 0, 5},  // 4

        {5, 0, 1}   // 5

    };

    // get the correct set of substrate vertices for this resolution

    CoordIJK *verts;

    if (isResolutionClassIII(*res))

        verts = vertsCIII;

    else

        verts = vertsCII;

    // adjust the center point to be in an aperture 33r substrate grid

    // these should be composed for speed

    _downAp3(&fijk->coord);

    _downAp3r(&fijk->coord);

    // if res is Class III we need to add a cw aperture 7 to get to

    // icosahedral Class II

    if (isResolutionClassIII(*res)) {

        _downAp7r(&fijk->coord);

        *res += 1;

    }

    // The center point is now in the same substrate grid as the origin

    // cell vertices. Add the center point substate coordinates

    // to each vertex to translate the vertices to that cell.

    for (int v = 0; v < NUM_HEX_VERTS; v++) {

        fijkVerts[v].face = fijk->face;

        _ijkAdd(&fijk->coord, &verts[v], &fijkVerts[v].coord);

        _ijkNormalize(&fijkVerts[v].coord);

    }

}

/**

 * Adjusts a FaceIJK address in place so that the resulting cell address is

 * relative to the correct icosahedral face.

 *

 * @param fijk The FaceIJK address of the cell.

 * @param res The H3 resolution of the cell.

 * @param pentLeading4 Whether or not the cell is a pentagon with a leading

 *        digit 4.

 * @param substrate Whether or not the cell is in a substrate grid.

 * @return 0 if on original face (no overage); 1 if on face edge (only occurs

 *         on substrate grids); 2 if overage on new face interior

 */

Overage _adjustOverageClassII(FaceIJK *fijk, int res, int pentLeading4,

                              int substrate) {

    Overage overage = NO_OVERAGE;

    CoordIJK *ijk = &fijk->coord;

    // get the maximum dimension value; scale if a substrate grid

    int maxDim = maxDimByCIIres[res];

    if (substrate) maxDim *= 3;

    // check for overage

    if (substrate && ijk->i + ijk->j + ijk->k == maxDim)  // on edge

        overage = FACE_EDGE;

    else if (ijk->i + ijk->j + ijk->k > maxDim)  // overage

    {

        overage = NEW_FACE;

        const FaceOrientIJK *fijkOrient;

        if (ijk->k > 0) {

            if (ijk->j > 0)  // jk "quadrant"

                fijkOrient = &faceNeighbors[fijk->face][JK];

            else  // ik "quadrant"

            {

                fijkOrient = &faceNeighbors[fijk->face][KI];

                // adjust for the pentagonal missing sequence

                if (pentLeading4) {

                    // translate origin to center of pentagon

                    CoordIJK origin;

                    _setIJK(&origin, maxDim, 0, 0);

                    CoordIJK tmp;

                    _ijkSub(ijk, &origin, &tmp);

                    // rotate to adjust for the missing sequence

                    _ijkRotate60cw(&tmp);

                    // translate the origin back to the center of the triangle

                    _ijkAdd(&tmp, &origin, ijk);

                }

            }

        } else  // ij "quadrant"

            fijkOrient = &faceNeighbors[fijk->face][IJ];

        fijk->face = fijkOrient->face;

        // rotate and translate for adjacent face

        for (int i = 0; i < fijkOrient->ccwRot60; i++) _ijkRotate60ccw(ijk);

        CoordIJK transVec = fijkOrient->translate;

        int unitScale = unitScaleByCIIres[res];

        if (substrate) unitScale *= 3;

        _ijkScale(&transVec, unitScale);

        _ijkAdd(ijk, &transVec, ijk);

        _ijkNormalize(ijk);

        // overage points on pentagon boundaries can end up on edges

        if (substrate && ijk->i + ijk->j + ijk->k == maxDim)  // on edge

            overage = FACE_EDGE;

    }

    return overage;

}

/**

 * Adjusts a FaceIJK address for a pentagon vertex in a substrate grid in

 * place so that the resulting cell address is relative to the correct

 * icosahedral face.

 *

 * @param fijk The FaceIJK address of the cell.

 * @param res The H3 resolution of the cell.

 */

Overage _adjustPentVertOverage(FaceIJK *fijk, int res) {

    int pentLeading4 = 0;

    Overage overage;

    do {

        overage = _adjustOverageClassII(fijk, res, pentLeading4, 1);

    } while (overage == NEW_FACE);

    return overage;

}

/**

 * Encodes a coordinate on the sphere to the corresponding icosahedral face and

 * containing the squared euclidean distance to that face center.

 *

 * Vec3d v is expected to be on the unit sphere.

 *

 * @param v The Vec3d coordinates to encode.

 * @param face Output: The icosahedral face containing the coordinates.

 * @param sqd Output: The squared euclidean distance to its face center.

 */

static void _vec3ToClosestFace(const Vec3d *v, int *face, double *sqd) {

    *face = 0;

    // The distance between two farthest points is 2.0, therefore the square of

    // the distance between two points should always be less or equal than 4.0 .

    *sqd = 5.0;

    for (int f = 0; f < NUM_ICOSA_FACES; ++f) {

        double sqdT = vec3DistSq(faceCenterPoint[f], *v);

        if (sqdT < *sqd) {

            *face = f;

            *sqd = sqdT;

        }

    }

}