-
-
Save user01/514a9ad748858310783c297d42659c45 to your computer and use it in GitHub Desktop.
pyopencl h3 example
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Number of platforms 1 | |
Platform Name NVIDIA CUDA | |
Platform Vendor NVIDIA Corporation | |
Platform Version OpenCL 3.0 CUDA 11.4.158 | |
Platform Profile FULL_PROFILE | |
Platform Extensions cl_khr_global_int32_base_atomics cl_khr_global_int32_extended_atomics cl_khr_local_int32_base_atomics cl_khr_local_int32_extended_atomics cl_khr_fp64 cl_khr_3d_image_writes cl_khr_byte_addressable_store cl_khr_icd cl_khr_gl_sharing cl_nv_compiler_options cl_nv_device_attribute_query cl_nv_pragma_unroll cl_nv_copy_opts cl_nv_create_buffer cl_khr_int64_base_atomics cl_khr_int64_extended_atomics cl_khr_device_uuid cl_khr_pci_bus_info | |
Platform Host timer resolution 0ns | |
Platform Extensions function suffix NV | |
Platform Name NVIDIA CUDA | |
Number of devices 1 | |
Device Name NVIDIA GeForce RTX 3090 | |
Device Vendor NVIDIA Corporation | |
Device Vendor ID 0x10de | |
Device Version OpenCL 3.0 CUDA | |
Driver Version 470.86 | |
Device OpenCL C Version OpenCL C 1.2 | |
Device Type GPU | |
Device Topology (NV) PCI-E, 01:00.0 | |
Device Profile FULL_PROFILE | |
Device Available Yes | |
Compiler Available Yes | |
Linker Available Yes | |
Max compute units 82 | |
Max clock frequency 1695MHz | |
Compute Capability (NV) 8.6 | |
Device Partition (core) | |
Max number of sub-devices 1 | |
Supported partition types None | |
Supported affinity domains (n/a) | |
Max work item dimensions 3 | |
Max work item sizes 1024x1024x64 | |
Max work group size 1024 | |
Preferred work group size multiple 32 | |
Warp size (NV) 32 | |
Max sub-groups per work group 0 | |
Preferred / native vector sizes | |
char 1 / 1 | |
short 1 / 1 | |
int 1 / 1 | |
long 1 / 1 | |
half 0 / 0 (n/a) | |
float 1 / 1 | |
double 1 / 1 (cl_khr_fp64) | |
Half-precision Floating-point support (n/a) | |
Single-precision Floating-point support (core) | |
Denormals Yes | |
Infinity and NANs Yes | |
Round to nearest Yes | |
Round to zero Yes | |
Round to infinity Yes | |
IEEE754-2008 fused multiply-add Yes | |
Support is emulated in software No | |
Correctly-rounded divide and sqrt operations Yes | |
Double-precision Floating-point support (cl_khr_fp64) | |
Denormals Yes | |
Infinity and NANs Yes | |
Round to nearest Yes | |
Round to zero Yes | |
Round to infinity Yes | |
IEEE754-2008 fused multiply-add Yes | |
Support is emulated in software No | |
Address bits 64, Little-Endian | |
Global memory size 25411846144 (23.67GiB) | |
Error Correction support No | |
Max memory allocation 6352961536 (5.917GiB) | |
Unified memory for Host and Device No | |
Integrated memory (NV) No | |
Shared Virtual Memory (SVM) capabilities (core) | |
Coarse-grained buffer sharing Yes | |
Fine-grained buffer sharing No | |
Fine-grained system sharing No | |
Atomics No | |
Minimum alignment for any data type 128 bytes | |
Alignment of base address 4096 bits (512 bytes) | |
Preferred alignment for atomics | |
SVM 0 bytes | |
Global 0 bytes | |
Local 0 bytes | |
Max size for global variable 0 | |
Preferred total size of global vars 0 | |
Global Memory cache type Read/Write | |
Global Memory cache size 2351104 (2.242MiB) | |
Global Memory cache line size 128 bytes | |
Image support Yes | |
Max number of samplers per kernel 32 | |
Max size for 1D images from buffer 268435456 pixels | |
Max 1D or 2D image array size 2048 images | |
Max 2D image size 32768x32768 pixels | |
Max 3D image size 16384x16384x16384 pixels | |
Max number of read image args 256 | |
Max number of write image args 32 | |
Max number of read/write image args 0 | |
Max number of pipe args 0 | |
Max active pipe reservations 0 | |
Max pipe packet size 0 | |
Local memory type Local | |
Local memory size 49152 (48KiB) | |
Registers per block (NV) 65536 | |
Max number of constant args 9 | |
Max constant buffer size 65536 (64KiB) | |
Max size of kernel argument 4352 (4.25KiB) | |
Queue properties (on host) | |
Out-of-order execution Yes | |
Profiling Yes | |
Queue properties (on device) | |
Out-of-order execution No | |
Profiling No | |
Preferred size 0 | |
Max size 0 | |
Max queues on device 0 | |
Max events on device 0 | |
Prefer user sync for interop No | |
Profiling timer resolution 1000ns | |
Execution capabilities | |
Run OpenCL kernels Yes | |
Run native kernels No | |
Sub-group independent forward progress No | |
Kernel execution timeout (NV) Yes | |
Concurrent copy and kernel execution (NV) Yes | |
Number of async copy engines 2 | |
IL version (n/a) | |
printf() buffer size 1048576 (1024KiB) | |
Built-in kernels (n/a) | |
Device Extensions cl_khr_global_int32_base_atomics cl_khr_global_int32_extended_atomics cl_khr_local_int32_base_atomics cl_khr_local_int32_extended_atomics cl_khr_fp64 cl_khr_3d_image_writes cl_khr_byte_addressable_store cl_khr_icd cl_khr_gl_sharing cl_nv_compiler_options cl_nv_device_attribute_query cl_nv_pragma_unroll cl_nv_copy_opts cl_nv_create_buffer cl_khr_int64_base_atomics cl_khr_int64_extended_atomics cl_khr_device_uuid cl_khr_pci_bus_info | |
NULL platform behavior | |
clGetPlatformInfo(NULL, CL_PLATFORM_NAME, ...) No platform | |
clGetDeviceIDs(NULL, CL_DEVICE_TYPE_ALL, ...) No platform | |
clCreateContext(NULL, ...) [default] No platform | |
clCreateContext(NULL, ...) [other] Success [NV] | |
clCreateContextFromType(NULL, CL_DEVICE_TYPE_DEFAULT) No platform | |
clCreateContextFromType(NULL, CL_DEVICE_TYPE_CPU) No devices found in platform | |
clCreateContextFromType(NULL, CL_DEVICE_TYPE_GPU) No platform | |
clCreateContextFromType(NULL, CL_DEVICE_TYPE_ACCELERATOR) No devices found in platform | |
clCreateContextFromType(NULL, CL_DEVICE_TYPE_CUSTOM) Invalid device type for platform | |
clCreateContextFromType(NULL, CL_DEVICE_TYPE_ALL) No platform | |
NOTE: your OpenCL library only supports OpenCL 2.0, | |
but some installed platforms support OpenCL 3.0. | |
Programs using 3.0 features may crash | |
or behave unexpectedly |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** Adapted from on https://github.com/uber/h3/tree/ee292a96906767162f2f530b46fc2428b2555748 */ | |
/** License Apache License Version 2.0, January 2004 */ | |
/** https://github.com/uber/h3/blob/ee292a96906767162f2f530b46fc2428b2555748/LICENSE */ | |
#define UINT64_C(value) ((ulong)(value)) | |
/** pi / 180 */ | |
#define M_PI_180 0.0174532925199432957692369076848861271111 | |
/** max H3 resolution; H3 version 1 has 16 resolutions, numbered 0 through 15 */ | |
#define MAX_H3_RES 15 | |
/** The number of faces on an icosahedron */ | |
#define NUM_ICOSA_FACES 20 | |
/** threshold epsilon */ | |
#define EPSILON 0.0000000000000001 | |
/** scaling factor from hex2d resolution 0 unit length | |
* (or distance between adjacent cell center points | |
* on the plane) to gnomonic unit length. */ | |
#define RES0_U_GNOMONIC 0.38196601125010500003 | |
/** 2.0 * PI */ | |
#define M_2PI 6.28318530717958647692528676655900576839433 | |
/** rotation angle between Class II and Class III resolution axes | |
* (asin(sqrt(3.0 / 28.0))) */ | |
#define M_AP7_ROT_RADS 0.333473172251832115336090755351601070065900389 | |
/** square root of 7 */ | |
#define M_SQRT7 2.6457513110645905905016157536392604257102 | |
/** sqrt(3) / 2.0 */ | |
#define M_SQRT3_2 0.8660254037844386467637231707529361834714 | |
/** sin(60') */ | |
#define M_SIN60 M_SQRT3_2 | |
/** The bit offset of the mode in an H3 index. */ | |
#define H3_MODE_OFFSET 59 | |
/** 1's in the 4 mode bits, 0's everywhere else. */ | |
#define H3_MODE_MASK ((ulong)(15) << H3_MODE_OFFSET) | |
/** 0's in the 4 mode bits, 1's everywhere else. */ | |
#define H3_MODE_MASK_NEGATIVE (~H3_MODE_MASK) | |
/** The bit offset of the resolution in an H3 index. */ | |
#define H3_RES_OFFSET 52 | |
/** 1's in the 4 resolution bits, 0's everywhere else. */ | |
#define H3_RES_MASK (UINT64_C(15) << H3_RES_OFFSET) | |
/** 0's in the 4 resolution bits, 1's everywhere else. */ | |
#define H3_RES_MASK_NEGATIVE (~H3_RES_MASK) | |
/** The bit offset of the base cell in an H3 index. */ | |
#define H3_BC_OFFSET 45 | |
/** 1's in the 7 base cell bits, 0's everywhere else. */ | |
#define H3_BC_MASK ((ulong)(127) << H3_BC_OFFSET) | |
/** 0's in the 7 base cell bits, 1's everywhere else. */ | |
#define H3_BC_MASK_NEGATIVE (~H3_BC_MASK) | |
#define H3_CELL_MODE 1 | |
/** Maximum input for any component to face-to-base-cell lookup functions */ | |
#define MAX_FACE_COORD 2 | |
/** | |
* Invalid index used to indicate an error from latLngToCell and related | |
* functions or missing data in arrays of H3 indices. Analogous to NaN in | |
* floating point. | |
*/ | |
#define H3_NULL 0 | |
/** The number of H3 base cells */ | |
#define NUM_BASE_CELLS 122 | |
/** | |
* Gets the integer resolution of h3. | |
*/ | |
#define H3_GET_RESOLUTION(h3) ((int)((((h3)&H3_RES_MASK) >> H3_RES_OFFSET))) | |
/** | |
* Sets the integer mode of h3 to v. | |
*/ | |
#define H3_SET_MODE(h3, v) \ | |
(h3) = (((h3)&H3_MODE_MASK_NEGATIVE) | (((ulong)(v)) << H3_MODE_OFFSET)) | |
/** | |
* Sets the integer resolution of h3. | |
*/ | |
#define H3_SET_RESOLUTION(h3, res) \ | |
(h3) = (((h3)&H3_RES_MASK_NEGATIVE) | (((ulong)(res)) << H3_RES_OFFSET)) | |
/** | |
* Sets the integer base cell of h3 to bc. | |
*/ | |
#define H3_SET_BASE_CELL(h3, bc) \ | |
(h3) = (((h3)&H3_BC_MASK_NEGATIVE) | (((ulong)(bc)) << H3_BC_OFFSET)) | |
/** | |
* H3 index with mode 0, res 0, base cell 0, and 7 for all index digits. | |
* Typically used to initialize the creation of an H3 cell index, which | |
* expects all direction digits to be 7 beyond the cell's resolution. | |
*/ | |
#define H3_INIT (UINT64_C(35184372088831)) | |
/** The number of bits in a single H3 resolution digit. */ | |
#define H3_PER_DIGIT_OFFSET 3 | |
/** 1's in the 3 bits of res 15 digit bits, 0's everywhere else. */ | |
#define H3_DIGIT_MASK ((ulong)(7)) | |
/** | |
* Sets the resolution res digit of h3 to the integer digit (0-7) | |
*/ | |
#define H3_SET_INDEX_DIGIT(h3, res, digit) \ | |
(h3) = \ | |
(((h3) & \ | |
~((H3_DIGIT_MASK << ((MAX_H3_RES - (res)) * H3_PER_DIGIT_OFFSET)))) | \ | |
(((ulong)(digit)) << ((MAX_H3_RES - (res)) * H3_PER_DIGIT_OFFSET))) | |
/** @brief H3 digit representing ijk+ axes direction. | |
* Values will be within the lowest 3 bits of an integer. | |
*/ | |
typedef enum { | |
/** H3 digit in center */ | |
CENTER_DIGIT = 0, | |
/** H3 digit in k-axes direction */ | |
K_AXES_DIGIT = 1, | |
/** H3 digit in j-axes direction */ | |
J_AXES_DIGIT = 2, | |
/** H3 digit in j == k direction */ | |
JK_AXES_DIGIT = J_AXES_DIGIT | K_AXES_DIGIT, /* 3 */ | |
/** H3 digit in i-axes direction */ | |
I_AXES_DIGIT = 4, | |
/** H3 digit in i == k direction */ | |
IK_AXES_DIGIT = I_AXES_DIGIT | K_AXES_DIGIT, /* 5 */ | |
/** H3 digit in i == j direction */ | |
IJ_AXES_DIGIT = I_AXES_DIGIT | J_AXES_DIGIT, /* 6 */ | |
/** H3 digit in the invalid direction */ | |
INVALID_DIGIT = 7, | |
/** Valid digits will be less than this value. Same value as INVALID_DIGIT. | |
*/ | |
NUM_DIGITS = INVALID_DIGIT, | |
/** Child digit which is skipped for pentagons */ | |
PENTAGON_SKIPPED_DIGIT = K_AXES_DIGIT /* 1 */ | |
} Direction; | |
/** @struct LatLng | |
@brief latitude/longitude in radians | |
*/ | |
typedef struct { | |
double lat; ///< latitude in radians | |
double lng; ///< longitude in radians | |
} LatLng; | |
/** @brief Result code (success or specific error) from an H3 operation */ | |
typedef uint H3Error; | |
/** @brief Identifier for an object (cell, edge, etc) in the H3 system. | |
* | |
* The H3Index fits within a 64-bit unsigned integer. | |
*/ | |
typedef ulong H3Index; | |
typedef enum { | |
E_SUCCESS = 0, // Success (no error) | |
E_FAILED = | |
1, // The operation failed but a more specific error is not available | |
E_DOMAIN = 2, // Argument was outside of acceptable range (when a more | |
// specific error code is not available) | |
E_LATLNG_DOMAIN = | |
3, // Latitude or longitude arguments were outside of acceptable range | |
E_RES_DOMAIN = 4, // Resolution argument was outside of acceptable range | |
E_CELL_INVALID = 5, // `H3Index` cell argument was not valid | |
E_DIR_EDGE_INVALID = 6, // `H3Index` directed edge argument was not valid | |
E_UNDIR_EDGE_INVALID = 7, // `H3Index` undirected edge argument was not valid | |
E_VERTEX_INVALID = 8, // `H3Index` vertex argument was not valid | |
E_PENTAGON = 9, // Pentagon distortion was encountered which the algorithm | |
// could not handle it | |
E_DUPLICATE_INPUT = 10, // Duplicate input was encountered in the arguments | |
// and the algorithm could not handle it | |
E_NOT_NEIGHBORS = 11, // `H3Index` cell arguments were not neighbors | |
E_RES_MISMATCH = 12, // `H3Index` cell arguments had incompatible resolutions | |
E_MEMORY = 13, // Necessary memory allocation failed | |
E_MEMORY_BOUNDS = 14, // Bounds of provided memory were not large enough | |
E_OPTION_INVALID = 15 // Mode or flags argument was not valid. | |
} H3ErrorCodes; | |
/** @struct CoordIJK | |
* @brief IJK hexagon coordinates | |
* | |
* Each axis is spaced 120 degrees apart. | |
*/ | |
typedef struct { | |
int i; ///< i component | |
int j; ///< j component | |
int k; ///< k component | |
} CoordIJK; | |
/** @struct FaceIJK | |
* @brief Face number and ijk coordinates on that face-centered coordinate | |
* system | |
*/ | |
typedef struct { | |
int face; ///< face number | |
CoordIJK coord; ///< ijk coordinates on that face | |
} FaceIJK; | |
/** @struct Vec2d | |
* @brief 2D floating-point vector | |
*/ | |
typedef struct { | |
double x; ///< x component | |
double y; ///< y component | |
} Vec2d; | |
/** @struct Vec3D | |
* @brief 3D floating point structure | |
*/ | |
typedef struct { | |
double x; ///< x component | |
double y; ///< y component | |
double z; ///< z component | |
} Vec3d; | |
/** @struct BaseCellRotation | |
* @brief base cell at a given ijk and required rotations into its system | |
*/ | |
typedef struct { | |
int baseCell; ///< base cell number | |
int ccwRot60; ///< number of ccw 60 degree rotations relative to current | |
/// face | |
} BaseCellRotation; | |
/** @struct BaseCellData | |
* @brief information on a single base cell | |
*/ | |
typedef struct { | |
FaceIJK homeFijk; ///< "home" face and normalized ijk coordinates on that face | |
int isPentagon; ///< is this base cell a pentagon? | |
int cwOffsetPent[2]; ///< if a pentagon, what are its two clockwise offset | |
/// faces? | |
} BaseCellData; | |
/** | |
* Gets the resolution res integer digit (0-7) of h3. | |
*/ | |
#define H3_GET_INDEX_DIGIT(h3, res) \ | |
((Direction)((((h3) >> ((MAX_H3_RES - (res)) * H3_PER_DIGIT_OFFSET)) & \ | |
H3_DIGIT_MASK))) | |
/** @brief CoordIJK unit vectors corresponding to the 7 H3 digits. | |
*/ | |
static const CoordIJK UNIT_VECS[] = { | |
{0, 0, 0}, // direction 0 | |
{0, 0, 1}, // direction 1 | |
{0, 1, 0}, // direction 2 | |
{0, 1, 1}, // direction 3 | |
{1, 0, 0}, // direction 4 | |
{1, 0, 1}, // direction 5 | |
{1, 1, 0} // direction 6 | |
}; | |
/** @brief Resolution 0 base cell data table. | |
* | |
* For each base cell, gives the "home" face and ijk+ coordinates on that face, | |
* whether or not the base cell is a pentagon. Additionally, if the base cell | |
* is a pentagon, the two cw offset rotation adjacent faces are given (-1 | |
* indicates that no cw offset rotation faces exist for this base cell). | |
*/ | |
const BaseCellData baseCellData[NUM_BASE_CELLS] = { | |
{{1, {1, 0, 0}}, 0, {0, 0}}, // base cell 0 | |
{{2, {1, 1, 0}}, 0, {0, 0}}, // base cell 1 | |
{{1, {0, 0, 0}}, 0, {0, 0}}, // base cell 2 | |
{{2, {1, 0, 0}}, 0, {0, 0}}, // base cell 3 | |
{{0, {2, 0, 0}}, 1, {-1, -1}}, // base cell 4 | |
{{1, {1, 1, 0}}, 0, {0, 0}}, // base cell 5 | |
{{1, {0, 0, 1}}, 0, {0, 0}}, // base cell 6 | |
{{2, {0, 0, 0}}, 0, {0, 0}}, // base cell 7 | |
{{0, {1, 0, 0}}, 0, {0, 0}}, // base cell 8 | |
{{2, {0, 1, 0}}, 0, {0, 0}}, // base cell 9 | |
{{1, {0, 1, 0}}, 0, {0, 0}}, // base cell 10 | |
{{1, {0, 1, 1}}, 0, {0, 0}}, // base cell 11 | |
{{3, {1, 0, 0}}, 0, {0, 0}}, // base cell 12 | |
{{3, {1, 1, 0}}, 0, {0, 0}}, // base cell 13 | |
{{11, {2, 0, 0}}, 1, {2, 6}}, // base cell 14 | |
{{4, {1, 0, 0}}, 0, {0, 0}}, // base cell 15 | |
{{0, {0, 0, 0}}, 0, {0, 0}}, // base cell 16 | |
{{6, {0, 1, 0}}, 0, {0, 0}}, // base cell 17 | |
{{0, {0, 0, 1}}, 0, {0, 0}}, // base cell 18 | |
{{2, {0, 1, 1}}, 0, {0, 0}}, // base cell 19 | |
{{7, {0, 0, 1}}, 0, {0, 0}}, // base cell 20 | |
{{2, {0, 0, 1}}, 0, {0, 0}}, // base cell 21 | |
{{0, {1, 1, 0}}, 0, {0, 0}}, // base cell 22 | |
{{6, {0, 0, 1}}, 0, {0, 0}}, // base cell 23 | |
{{10, {2, 0, 0}}, 1, {1, 5}}, // base cell 24 | |
{{6, {0, 0, 0}}, 0, {0, 0}}, // base cell 25 | |
{{3, {0, 0, 0}}, 0, {0, 0}}, // base cell 26 | |
{{11, {1, 0, 0}}, 0, {0, 0}}, // base cell 27 | |
{{4, {1, 1, 0}}, 0, {0, 0}}, // base cell 28 | |
{{3, {0, 1, 0}}, 0, {0, 0}}, // base cell 29 | |
{{0, {0, 1, 1}}, 0, {0, 0}}, // base cell 30 | |
{{4, {0, 0, 0}}, 0, {0, 0}}, // base cell 31 | |
{{5, {0, 1, 0}}, 0, {0, 0}}, // base cell 32 | |
{{0, {0, 1, 0}}, 0, {0, 0}}, // base cell 33 | |
{{7, {0, 1, 0}}, 0, {0, 0}}, // base cell 34 | |
{{11, {1, 1, 0}}, 0, {0, 0}}, // base cell 35 | |
{{7, {0, 0, 0}}, 0, {0, 0}}, // base cell 36 | |
{{10, {1, 0, 0}}, 0, {0, 0}}, // base cell 37 | |
{{12, {2, 0, 0}}, 1, {3, 7}}, // base cell 38 | |
{{6, {1, 0, 1}}, 0, {0, 0}}, // base cell 39 | |
{{7, {1, 0, 1}}, 0, {0, 0}}, // base cell 40 | |
{{4, {0, 0, 1}}, 0, {0, 0}}, // base cell 41 | |
{{3, {0, 0, 1}}, 0, {0, 0}}, // base cell 42 | |
{{3, {0, 1, 1}}, 0, {0, 0}}, // base cell 43 | |
{{4, {0, 1, 0}}, 0, {0, 0}}, // base cell 44 | |
{{6, {1, 0, 0}}, 0, {0, 0}}, // base cell 45 | |
{{11, {0, 0, 0}}, 0, {0, 0}}, // base cell 46 | |
{{8, {0, 0, 1}}, 0, {0, 0}}, // base cell 47 | |
{{5, {0, 0, 1}}, 0, {0, 0}}, // base cell 48 | |
{{14, {2, 0, 0}}, 1, {0, 9}}, // base cell 49 | |
{{5, {0, 0, 0}}, 0, {0, 0}}, // base cell 50 | |
{{12, {1, 0, 0}}, 0, {0, 0}}, // base cell 51 | |
{{10, {1, 1, 0}}, 0, {0, 0}}, // base cell 52 | |
{{4, {0, 1, 1}}, 0, {0, 0}}, // base cell 53 | |
{{12, {1, 1, 0}}, 0, {0, 0}}, // base cell 54 | |
{{7, {1, 0, 0}}, 0, {0, 0}}, // base cell 55 | |
{{11, {0, 1, 0}}, 0, {0, 0}}, // base cell 56 | |
{{10, {0, 0, 0}}, 0, {0, 0}}, // base cell 57 | |
{{13, {2, 0, 0}}, 1, {4, 8}}, // base cell 58 | |
{{10, {0, 0, 1}}, 0, {0, 0}}, // base cell 59 | |
{{11, {0, 0, 1}}, 0, {0, 0}}, // base cell 60 | |
{{9, {0, 1, 0}}, 0, {0, 0}}, // base cell 61 | |
{{8, {0, 1, 0}}, 0, {0, 0}}, // base cell 62 | |
{{6, {2, 0, 0}}, 1, {11, 15}}, // base cell 63 | |
{{8, {0, 0, 0}}, 0, {0, 0}}, // base cell 64 | |
{{9, {0, 0, 1}}, 0, {0, 0}}, // base cell 65 | |
{{14, {1, 0, 0}}, 0, {0, 0}}, // base cell 66 | |
{{5, {1, 0, 1}}, 0, {0, 0}}, // base cell 67 | |
{{16, {0, 1, 1}}, 0, {0, 0}}, // base cell 68 | |
{{8, {1, 0, 1}}, 0, {0, 0}}, // base cell 69 | |
{{5, {1, 0, 0}}, 0, {0, 0}}, // base cell 70 | |
{{12, {0, 0, 0}}, 0, {0, 0}}, // base cell 71 | |
{{7, {2, 0, 0}}, 1, {12, 16}}, // base cell 72 | |
{{12, {0, 1, 0}}, 0, {0, 0}}, // base cell 73 | |
{{10, {0, 1, 0}}, 0, {0, 0}}, // base cell 74 | |
{{9, {0, 0, 0}}, 0, {0, 0}}, // base cell 75 | |
{{13, {1, 0, 0}}, 0, {0, 0}}, // base cell 76 | |
{{16, {0, 0, 1}}, 0, {0, 0}}, // base cell 77 | |
{{15, {0, 1, 1}}, 0, {0, 0}}, // base cell 78 | |
{{15, {0, 1, 0}}, 0, {0, 0}}, // base cell 79 | |
{{16, {0, 1, 0}}, 0, {0, 0}}, // base cell 80 | |
{{14, {1, 1, 0}}, 0, {0, 0}}, // base cell 81 | |
{{13, {1, 1, 0}}, 0, {0, 0}}, // base cell 82 | |
{{5, {2, 0, 0}}, 1, {10, 19}}, // base cell 83 | |
{{8, {1, 0, 0}}, 0, {0, 0}}, // base cell 84 | |
{{14, {0, 0, 0}}, 0, {0, 0}}, // base cell 85 | |
{{9, {1, 0, 1}}, 0, {0, 0}}, // base cell 86 | |
{{14, {0, 0, 1}}, 0, {0, 0}}, // base cell 87 | |
{{17, {0, 0, 1}}, 0, {0, 0}}, // base cell 88 | |
{{12, {0, 0, 1}}, 0, {0, 0}}, // base cell 89 | |
{{16, {0, 0, 0}}, 0, {0, 0}}, // base cell 90 | |
{{17, {0, 1, 1}}, 0, {0, 0}}, // base cell 91 | |
{{15, {0, 0, 1}}, 0, {0, 0}}, // base cell 92 | |
{{16, {1, 0, 1}}, 0, {0, 0}}, // base cell 93 | |
{{9, {1, 0, 0}}, 0, {0, 0}}, // base cell 94 | |
{{15, {0, 0, 0}}, 0, {0, 0}}, // base cell 95 | |
{{13, {0, 0, 0}}, 0, {0, 0}}, // base cell 96 | |
{{8, {2, 0, 0}}, 1, {13, 17}}, // base cell 97 | |
{{13, {0, 1, 0}}, 0, {0, 0}}, // base cell 98 | |
{{17, {1, 0, 1}}, 0, {0, 0}}, // base cell 99 | |
{{19, {0, 1, 0}}, 0, {0, 0}}, // base cell 100 | |
{{14, {0, 1, 0}}, 0, {0, 0}}, // base cell 101 | |
{{19, {0, 1, 1}}, 0, {0, 0}}, // base cell 102 | |
{{17, {0, 1, 0}}, 0, {0, 0}}, // base cell 103 | |
{{13, {0, 0, 1}}, 0, {0, 0}}, // base cell 104 | |
{{17, {0, 0, 0}}, 0, {0, 0}}, // base cell 105 | |
{{16, {1, 0, 0}}, 0, {0, 0}}, // base cell 106 | |
{{9, {2, 0, 0}}, 1, {14, 18}}, // base cell 107 | |
{{15, {1, 0, 1}}, 0, {0, 0}}, // base cell 108 | |
{{15, {1, 0, 0}}, 0, {0, 0}}, // base cell 109 | |
{{18, {0, 1, 1}}, 0, {0, 0}}, // base cell 110 | |
{{18, {0, 0, 1}}, 0, {0, 0}}, // base cell 111 | |
{{19, {0, 0, 1}}, 0, {0, 0}}, // base cell 112 | |
{{17, {1, 0, 0}}, 0, {0, 0}}, // base cell 113 | |
{{19, {0, 0, 0}}, 0, {0, 0}}, // base cell 114 | |
{{18, {0, 1, 0}}, 0, {0, 0}}, // base cell 115 | |
{{18, {1, 0, 1}}, 0, {0, 0}}, // base cell 116 | |
{{19, {2, 0, 0}}, 1, {-1, -1}}, // base cell 117 | |
{{19, {1, 0, 0}}, 0, {0, 0}}, // base cell 118 | |
{{18, {0, 0, 0}}, 0, {0, 0}}, // base cell 119 | |
{{19, {1, 0, 1}}, 0, {0, 0}}, // base cell 120 | |
{{18, {1, 0, 0}}, 0, {0, 0}} // base cell 121 | |
}; | |
/** @brief icosahedron face centers in lat/lng radians */ | |
const LatLng faceCenterGeo[NUM_ICOSA_FACES] = { | |
{0.803582649718989942, 1.248397419617396099}, // face 0 | |
{1.307747883455638156, 2.536945009877921159}, // face 1 | |
{1.054751253523952054, -1.347517358900396623}, // face 2 | |
{0.600191595538186799, -0.450603909469755746}, // face 3 | |
{0.491715428198773866, 0.401988202911306943}, // face 4 | |
{0.172745327415618701, 1.678146885280433686}, // face 5 | |
{0.605929321571350690, 2.953923329812411617}, // face 6 | |
{0.427370518328979641, -1.888876200336285401}, // face 7 | |
{-0.079066118549212831, -0.733429513380867741}, // face 8 | |
{-0.230961644455383637, 0.506495587332349035}, // face 9 | |
{0.079066118549212831, 2.408163140208925497}, // face 10 | |
{0.230961644455383637, -2.635097066257444203}, // face 11 | |
{-0.172745327415618701, -1.463445768309359553}, // face 12 | |
{-0.605929321571350690, -0.187669323777381622}, // face 13 | |
{-0.427370518328979641, 1.252716453253507838}, // face 14 | |
{-0.600191595538186799, 2.690988744120037492}, // face 15 | |
{-0.491715428198773866, -2.739604450678486295}, // face 16 | |
{-0.803582649718989942, -1.893195233972397139}, // face 17 | |
{-1.307747883455638156, -0.604647643711872080}, // face 18 | |
{-1.054751253523952054, 1.794075294689396615}, // face 19 | |
}; | |
/** @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 Resolution 0 base cell lookup table for each face. | |
* | |
* Given the face number and a resolution 0 ijk+ coordinate in that face's | |
* face-centered ijk coordinate system, gives the base cell located at that | |
* coordinate and the number of 60 ccw rotations to rotate into that base | |
* cell's orientation. | |
* | |
* Valid lookup coordinates are from (0, 0, 0) to (2, 2, 2). | |
* | |
* This table can be accessed using the functions `_faceIjkToBaseCell` and | |
* `_faceIjkToBaseCellCCWrot60` | |
*/ | |
static const BaseCellRotation | |
faceIjkBaseCells[NUM_ICOSA_FACES][3][3][3] = { | |
{// face 0 | |
{ | |
// i 0 | |
{{16, 0}, {18, 0}, {24, 0}}, // j 0 | |
{{33, 0}, {30, 0}, {32, 3}}, // j 1 | |
{{49, 1}, {48, 3}, {50, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{8, 0}, {5, 5}, {10, 5}}, // j 0 | |
{{22, 0}, {16, 0}, {18, 0}}, // j 1 | |
{{41, 1}, {33, 0}, {30, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{4, 0}, {0, 5}, {2, 5}}, // j 0 | |
{{15, 1}, {8, 0}, {5, 5}}, // j 1 | |
{{31, 1}, {22, 0}, {16, 0}} // j 2 | |
}}, | |
{// face 1 | |
{ | |
// i 0 | |
{{2, 0}, {6, 0}, {14, 0}}, // j 0 | |
{{10, 0}, {11, 0}, {17, 3}}, // j 1 | |
{{24, 1}, {23, 3}, {25, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{0, 0}, {1, 5}, {9, 5}}, // j 0 | |
{{5, 0}, {2, 0}, {6, 0}}, // j 1 | |
{{18, 1}, {10, 0}, {11, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{4, 1}, {3, 5}, {7, 5}}, // j 0 | |
{{8, 1}, {0, 0}, {1, 5}}, // j 1 | |
{{16, 1}, {5, 0}, {2, 0}} // j 2 | |
}}, | |
{// face 2 | |
{ | |
// i 0 | |
{{7, 0}, {21, 0}, {38, 0}}, // j 0 | |
{{9, 0}, {19, 0}, {34, 3}}, // j 1 | |
{{14, 1}, {20, 3}, {36, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{3, 0}, {13, 5}, {29, 5}}, // j 0 | |
{{1, 0}, {7, 0}, {21, 0}}, // j 1 | |
{{6, 1}, {9, 0}, {19, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{4, 2}, {12, 5}, {26, 5}}, // j 0 | |
{{0, 1}, {3, 0}, {13, 5}}, // j 1 | |
{{2, 1}, {1, 0}, {7, 0}} // j 2 | |
}}, | |
{// face 3 | |
{ | |
// i 0 | |
{{26, 0}, {42, 0}, {58, 0}}, // j 0 | |
{{29, 0}, {43, 0}, {62, 3}}, // j 1 | |
{{38, 1}, {47, 3}, {64, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{12, 0}, {28, 5}, {44, 5}}, // j 0 | |
{{13, 0}, {26, 0}, {42, 0}}, // j 1 | |
{{21, 1}, {29, 0}, {43, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{4, 3}, {15, 5}, {31, 5}}, // j 0 | |
{{3, 1}, {12, 0}, {28, 5}}, // j 1 | |
{{7, 1}, {13, 0}, {26, 0}} // j 2 | |
}}, | |
{// face 4 | |
{ | |
// i 0 | |
{{31, 0}, {41, 0}, {49, 0}}, // j 0 | |
{{44, 0}, {53, 0}, {61, 3}}, // j 1 | |
{{58, 1}, {65, 3}, {75, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{15, 0}, {22, 5}, {33, 5}}, // j 0 | |
{{28, 0}, {31, 0}, {41, 0}}, // j 1 | |
{{42, 1}, {44, 0}, {53, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{4, 4}, {8, 5}, {16, 5}}, // j 0 | |
{{12, 1}, {15, 0}, {22, 5}}, // j 1 | |
{{26, 1}, {28, 0}, {31, 0}} // j 2 | |
}}, | |
{// face 5 | |
{ | |
// i 0 | |
{{50, 0}, {48, 0}, {49, 3}}, // j 0 | |
{{32, 0}, {30, 3}, {33, 3}}, // j 1 | |
{{24, 3}, {18, 3}, {16, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{70, 0}, {67, 0}, {66, 3}}, // j 0 | |
{{52, 3}, {50, 0}, {48, 0}}, // j 1 | |
{{37, 3}, {32, 0}, {30, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{83, 0}, {87, 3}, {85, 3}}, // j 0 | |
{{74, 3}, {70, 0}, {67, 0}}, // j 1 | |
{{57, 1}, {52, 3}, {50, 0}} // j 2 | |
}}, | |
{// face 6 | |
{ | |
// i 0 | |
{{25, 0}, {23, 0}, {24, 3}}, // j 0 | |
{{17, 0}, {11, 3}, {10, 3}}, // j 1 | |
{{14, 3}, {6, 3}, {2, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{45, 0}, {39, 0}, {37, 3}}, // j 0 | |
{{35, 3}, {25, 0}, {23, 0}}, // j 1 | |
{{27, 3}, {17, 0}, {11, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{63, 0}, {59, 3}, {57, 3}}, // j 0 | |
{{56, 3}, {45, 0}, {39, 0}}, // j 1 | |
{{46, 3}, {35, 3}, {25, 0}} // j 2 | |
}}, | |
{// face 7 | |
{ | |
// i 0 | |
{{36, 0}, {20, 0}, {14, 3}}, // j 0 | |
{{34, 0}, {19, 3}, {9, 3}}, // j 1 | |
{{38, 3}, {21, 3}, {7, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{55, 0}, {40, 0}, {27, 3}}, // j 0 | |
{{54, 3}, {36, 0}, {20, 0}}, // j 1 | |
{{51, 3}, {34, 0}, {19, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{72, 0}, {60, 3}, {46, 3}}, // j 0 | |
{{73, 3}, {55, 0}, {40, 0}}, // j 1 | |
{{71, 3}, {54, 3}, {36, 0}} // j 2 | |
}}, | |
{// face 8 | |
{ | |
// i 0 | |
{{64, 0}, {47, 0}, {38, 3}}, // j 0 | |
{{62, 0}, {43, 3}, {29, 3}}, // j 1 | |
{{58, 3}, {42, 3}, {26, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{84, 0}, {69, 0}, {51, 3}}, // j 0 | |
{{82, 3}, {64, 0}, {47, 0}}, // j 1 | |
{{76, 3}, {62, 0}, {43, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{97, 0}, {89, 3}, {71, 3}}, // j 0 | |
{{98, 3}, {84, 0}, {69, 0}}, // j 1 | |
{{96, 3}, {82, 3}, {64, 0}} // j 2 | |
}}, | |
{// face 9 | |
{ | |
// i 0 | |
{{75, 0}, {65, 0}, {58, 3}}, // j 0 | |
{{61, 0}, {53, 3}, {44, 3}}, // j 1 | |
{{49, 3}, {41, 3}, {31, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{94, 0}, {86, 0}, {76, 3}}, // j 0 | |
{{81, 3}, {75, 0}, {65, 0}}, // j 1 | |
{{66, 3}, {61, 0}, {53, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{107, 0}, {104, 3}, {96, 3}}, // j 0 | |
{{101, 3}, {94, 0}, {86, 0}}, // j 1 | |
{{85, 3}, {81, 3}, {75, 0}} // j 2 | |
}}, | |
{// face 10 | |
{ | |
// i 0 | |
{{57, 0}, {59, 0}, {63, 3}}, // j 0 | |
{{74, 0}, {78, 3}, {79, 3}}, // j 1 | |
{{83, 3}, {92, 3}, {95, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{37, 0}, {39, 3}, {45, 3}}, // j 0 | |
{{52, 0}, {57, 0}, {59, 0}}, // j 1 | |
{{70, 3}, {74, 0}, {78, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{24, 0}, {23, 3}, {25, 3}}, // j 0 | |
{{32, 3}, {37, 0}, {39, 3}}, // j 1 | |
{{50, 3}, {52, 0}, {57, 0}} // j 2 | |
}}, | |
{// face 11 | |
{ | |
// i 0 | |
{{46, 0}, {60, 0}, {72, 3}}, // j 0 | |
{{56, 0}, {68, 3}, {80, 3}}, // j 1 | |
{{63, 3}, {77, 3}, {90, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{27, 0}, {40, 3}, {55, 3}}, // j 0 | |
{{35, 0}, {46, 0}, {60, 0}}, // j 1 | |
{{45, 3}, {56, 0}, {68, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{14, 0}, {20, 3}, {36, 3}}, // j 0 | |
{{17, 3}, {27, 0}, {40, 3}}, // j 1 | |
{{25, 3}, {35, 0}, {46, 0}} // j 2 | |
}}, | |
{// face 12 | |
{ | |
// i 0 | |
{{71, 0}, {89, 0}, {97, 3}}, // j 0 | |
{{73, 0}, {91, 3}, {103, 3}}, // j 1 | |
{{72, 3}, {88, 3}, {105, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{51, 0}, {69, 3}, {84, 3}}, // j 0 | |
{{54, 0}, {71, 0}, {89, 0}}, // j 1 | |
{{55, 3}, {73, 0}, {91, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{38, 0}, {47, 3}, {64, 3}}, // j 0 | |
{{34, 3}, {51, 0}, {69, 3}}, // j 1 | |
{{36, 3}, {54, 0}, {71, 0}} // j 2 | |
}}, | |
{// face 13 | |
{ | |
// i 0 | |
{{96, 0}, {104, 0}, {107, 3}}, // j 0 | |
{{98, 0}, {110, 3}, {115, 3}}, // j 1 | |
{{97, 3}, {111, 3}, {119, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{76, 0}, {86, 3}, {94, 3}}, // j 0 | |
{{82, 0}, {96, 0}, {104, 0}}, // j 1 | |
{{84, 3}, {98, 0}, {110, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{58, 0}, {65, 3}, {75, 3}}, // j 0 | |
{{62, 3}, {76, 0}, {86, 3}}, // j 1 | |
{{64, 3}, {82, 0}, {96, 0}} // j 2 | |
}}, | |
{// face 14 | |
{ | |
// i 0 | |
{{85, 0}, {87, 0}, {83, 3}}, // j 0 | |
{{101, 0}, {102, 3}, {100, 3}}, // j 1 | |
{{107, 3}, {112, 3}, {114, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{66, 0}, {67, 3}, {70, 3}}, // j 0 | |
{{81, 0}, {85, 0}, {87, 0}}, // j 1 | |
{{94, 3}, {101, 0}, {102, 3}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{49, 0}, {48, 3}, {50, 3}}, // j 0 | |
{{61, 3}, {66, 0}, {67, 3}}, // j 1 | |
{{75, 3}, {81, 0}, {85, 0}} // j 2 | |
}}, | |
{// face 15 | |
{ | |
// i 0 | |
{{95, 0}, {92, 0}, {83, 0}}, // j 0 | |
{{79, 0}, {78, 0}, {74, 3}}, // j 1 | |
{{63, 1}, {59, 3}, {57, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{109, 0}, {108, 0}, {100, 5}}, // j 0 | |
{{93, 1}, {95, 0}, {92, 0}}, // j 1 | |
{{77, 1}, {79, 0}, {78, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{117, 4}, {118, 5}, {114, 5}}, // j 0 | |
{{106, 1}, {109, 0}, {108, 0}}, // j 1 | |
{{90, 1}, {93, 1}, {95, 0}} // j 2 | |
}}, | |
{// face 16 | |
{ | |
// i 0 | |
{{90, 0}, {77, 0}, {63, 0}}, // j 0 | |
{{80, 0}, {68, 0}, {56, 3}}, // j 1 | |
{{72, 1}, {60, 3}, {46, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{106, 0}, {93, 0}, {79, 5}}, // j 0 | |
{{99, 1}, {90, 0}, {77, 0}}, // j 1 | |
{{88, 1}, {80, 0}, {68, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{117, 3}, {109, 5}, {95, 5}}, // j 0 | |
{{113, 1}, {106, 0}, {93, 0}}, // j 1 | |
{{105, 1}, {99, 1}, {90, 0}} // j 2 | |
}}, | |
{// face 17 | |
{ | |
// i 0 | |
{{105, 0}, {88, 0}, {72, 0}}, // j 0 | |
{{103, 0}, {91, 0}, {73, 3}}, // j 1 | |
{{97, 1}, {89, 3}, {71, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{113, 0}, {99, 0}, {80, 5}}, // j 0 | |
{{116, 1}, {105, 0}, {88, 0}}, // j 1 | |
{{111, 1}, {103, 0}, {91, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{117, 2}, {106, 5}, {90, 5}}, // j 0 | |
{{121, 1}, {113, 0}, {99, 0}}, // j 1 | |
{{119, 1}, {116, 1}, {105, 0}} // j 2 | |
}}, | |
{// face 18 | |
{ | |
// i 0 | |
{{119, 0}, {111, 0}, {97, 0}}, // j 0 | |
{{115, 0}, {110, 0}, {98, 3}}, // j 1 | |
{{107, 1}, {104, 3}, {96, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{121, 0}, {116, 0}, {103, 5}}, // j 0 | |
{{120, 1}, {119, 0}, {111, 0}}, // j 1 | |
{{112, 1}, {115, 0}, {110, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{117, 1}, {113, 5}, {105, 5}}, // j 0 | |
{{118, 1}, {121, 0}, {116, 0}}, // j 1 | |
{{114, 1}, {120, 1}, {119, 0}} // j 2 | |
}}, | |
{// face 19 | |
{ | |
// i 0 | |
{{114, 0}, {112, 0}, {107, 0}}, // j 0 | |
{{100, 0}, {102, 0}, {101, 3}}, // j 1 | |
{{83, 1}, {87, 3}, {85, 3}} // j 2 | |
}, | |
{ | |
// i 1 | |
{{118, 0}, {120, 0}, {115, 5}}, // j 0 | |
{{108, 1}, {114, 0}, {112, 0}}, // j 1 | |
{{92, 1}, {100, 0}, {102, 0}} // j 2 | |
}, | |
{ | |
// i 2 | |
{{117, 0}, {121, 5}, {119, 5}}, // j 0 | |
{{109, 1}, {118, 0}, {120, 0}}, // j 1 | |
{{95, 1}, {108, 1}, {114, 0}} // j 2 | |
}}}; | |
/** | |
* Convert from decimal degrees to radians. | |
* | |
* @param degrees The decimal degrees. | |
* @return The corresponding radians. | |
*/ | |
double degsToRads(double degrees) { return degrees * M_PI_180; } | |
/** | |
* Normalizes radians to a value between 0.0 and two PI. | |
* | |
* @param rads The input radians value. | |
* @return The normalized radians value. | |
*/ | |
double _posAngleRads(double rads) { | |
double tmp = ((rads < 0.0) ? rads + M_2PI : rads); | |
if (rads >= M_2PI) | |
tmp -= M_2PI; | |
return tmp; | |
} | |
/** | |
* Determines the azimuth to p2 from p1 in radians. | |
* | |
* @param p1 The first spherical coordinates. | |
* @param p2 The second spherical coordinates. | |
* @return The azimuth in radians from p1 to p2. | |
*/ | |
double _geoAzimuthRads(const LatLng *p1, const LatLng *p2) { | |
return atan2(cos(p2->lat) * sin(p2->lng - p1->lng), | |
cos(p1->lat) * sin(p2->lat) - | |
sin(p1->lat) * cos(p2->lat) * cos(p2->lng - p1->lng)); | |
} | |
/** | |
* Returns whether or not a resolution is a Class III grid. Note that odd | |
* resolutions are Class III and even resolutions are Class II. | |
* @param res The H3 resolution. | |
* @return 1 if the resolution is a Class III grid, and 0 if the resolution is | |
* a Class II grid. | |
*/ | |
int isResolutionClassIII(int res) { return res % 2; } | |
/** | |
* Calculate the 3D coordinate on unit sphere from the latitude and longitude. | |
* | |
* @param geo The latitude and longitude of the point. | |
* @param v The 3D coordinate of the point. | |
*/ | |
void _geoToVec3d(const LatLng *geo, Vec3d *v) { | |
double r = cos(geo->lat); | |
v->z = sin(geo->lat); | |
v->x = cos(geo->lng) * r; | |
v->y = sin(geo->lng) * r; | |
} | |
/** | |
* Square of a number | |
* | |
* @param x The input number. | |
* @return The square of the input number. | |
*/ | |
double _square(double x) { return x * x; } | |
/** | |
* Calculate the square of the distance between two 3D coordinates. | |
* | |
* @param v1 The first 3D coordinate. | |
* @param v2 The second 3D coordinate. | |
* @return The square of the distance between the given points. | |
*/ | |
double _pointSquareDist(const Vec3d *v1, const Vec3d *v2) { | |
return _square(v1->x - v2->x) + _square(v1->y - v2->y) + | |
_square(v1->z - v2->z); | |
} | |
/** | |
* Encodes a coordinate on the sphere to the corresponding icosahedral face and | |
* containing 2D hex coordinates relative to that face center. | |
* | |
* @param g The spherical coordinates to encode. | |
* @param res The desired H3 resolution for the encoding. | |
* @param face The icosahedral face containing the spherical coordinates. | |
* @param v The 2D hex coordinates of the cell containing the point. | |
*/ | |
void _geoToHex2d(const LatLng *g, int res, int *face, Vec2d *v) { | |
Vec3d v3d; | |
_geoToVec3d(g, &v3d); | |
// determine the icosahedron face | |
*face = 0; | |
double sqd = _pointSquareDist(&faceCenterPoint[0], &v3d); | |
for (int f = 1; f < NUM_ICOSA_FACES; f++) { | |
double sqdT = _pointSquareDist(&faceCenterPoint[f], &v3d); | |
if (sqdT < sqd) { | |
*face = f; | |
sqd = sqdT; | |
} | |
} | |
// cos(r) = 1 - 2 * sin^2(r/2) = 1 - 2 * (sqd / 4) = 1 - sqd/2 | |
double r = acos(1 - sqd / 2); | |
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(_geoAzimuthRads(&faceCenterGeo[*face], g))); | |
// 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 /= 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); | |
} | |
/** | |
* Normalizes ijk coordinates by setting the components to the smallest possible | |
* values. Works in place. | |
* | |
* @param c The ijk coordinates to normalize. | |
*/ | |
void _ijkNormalize(CoordIJK *c) { | |
// remove any negative values | |
if (c->i < 0) { | |
c->j -= c->i; | |
c->k -= c->i; | |
c->i = 0; | |
} | |
if (c->j < 0) { | |
c->i -= c->j; | |
c->k -= c->j; | |
c->j = 0; | |
} | |
if (c->k < 0) { | |
c->i -= c->k; | |
c->j -= c->k; | |
c->k = 0; | |
} | |
// remove the min value if needed | |
int min = c->i; | |
if (c->j < min) | |
min = c->j; | |
if (c->k < min) | |
min = c->k; | |
if (min > 0) { | |
c->i -= min; | |
c->j -= min; | |
c->k -= min; | |
} | |
} | |
/** | |
* Determine the containing hex in ijk+ coordinates for a 2D cartesian | |
* coordinate vector (from DGGRID). | |
* | |
* @param v The 2D cartesian coordinate vector. | |
* @param h The ijk+ coordinates of the containing hex. | |
*/ | |
void _hex2dToCoordIJK(const Vec2d *v, CoordIJK *h) { | |
double a1, a2; | |
double x1, x2; | |
int m1, m2; | |
double r1, r2; | |
// quantize into the ij system and then normalize | |
h->k = 0; | |
a1 = fabs(v->x); // These used to be fabsl calls, but undefined on device | |
a2 = fabs(v->y); | |
// first do a reverse conversion | |
x2 = a2 / M_SIN60; | |
x1 = a1 + x2 / 2.0; | |
// check if we have the center of a hex | |
m1 = x1; | |
m2 = x2; | |
// otherwise round correctly | |
r1 = x1 - m1; | |
r2 = x2 - m2; | |
if (r1 < 0.5) { | |
if (r1 < 1.0 / 3.0) { | |
if (r2 < (1.0 + r1) / 2.0) { | |
h->i = m1; | |
h->j = m2; | |
} else { | |
h->i = m1; | |
h->j = m2 + 1; | |
} | |
} else { | |
if (r2 < (1.0 - r1)) { | |
h->j = m2; | |
} else { | |
h->j = m2 + 1; | |
} | |
if ((1.0 - r1) <= r2 && r2 < (2.0 * r1)) { | |
h->i = m1 + 1; | |
} else { | |
h->i = m1; | |
} | |
} | |
} else { | |
if (r1 < 2.0 / 3.0) { | |
if (r2 < (1.0 - r1)) { | |
h->j = m2; | |
} else { | |
h->j = m2 + 1; | |
} | |
if ((2.0 * r1 - 1.0) < r2 && r2 < (1.0 - r1)) { | |
h->i = m1; | |
} else { | |
h->i = m1 + 1; | |
} | |
} else { | |
if (r2 < (r1 / 2.0)) { | |
h->i = m1 + 1; | |
h->j = m2; | |
} else { | |
h->i = m1 + 1; | |
h->j = m2 + 1; | |
} | |
} | |
} | |
// now fold across the axes if necessary | |
if (v->x < 0.0) { | |
if ((h->j % 2) == 0) // even | |
{ | |
ulong axisi = h->j / 2; | |
ulong diff = h->i - axisi; | |
h->i = h->i - 2.0 * diff; | |
} else { | |
ulong axisi = (h->j + 1) / 2; | |
ulong diff = h->i - axisi; | |
h->i = h->i - (2.0 * diff + 1); | |
} | |
} | |
if (v->y < 0.0) { | |
h->i = h->i - (2 * h->j + 1) / 2; | |
h->j = -1 * h->j; | |
} | |
_ijkNormalize(h); | |
} | |
/** @brief Find base cell given FaceIJK. | |
* | |
* Given the face number and a resolution 0 ijk+ coordinate in that face's | |
* face-centered ijk coordinate system, return the base cell located at that | |
* coordinate. | |
* | |
* Valid ijk+ lookup coordinates are from (0, 0, 0) to (2, 2, 2). | |
*/ | |
int _faceIjkToBaseCell(const FaceIJK *h) { | |
return faceIjkBaseCells[h->face][h->coord.i][h->coord.j][h->coord.k].baseCell; | |
} | |
/** | |
* Encodes a coordinate on the sphere to the FaceIJK address of the containing | |
* cell at the specified resolution. | |
* | |
* @param g The spherical coordinates to encode. | |
* @param res The desired H3 resolution for the encoding. | |
* @param h The FaceIJK address of the containing cell at resolution res. | |
*/ | |
void _geoToFaceIjk(const LatLng *g, int res, FaceIJK *h) { | |
// first convert to hex2d | |
Vec2d v; | |
_geoToHex2d(g, res, &h->face, &v); | |
// then convert to ijk+ | |
_hex2dToCoordIJK(&v, &h->coord); | |
} | |
/** | |
* Add two ijk coordinates. | |
* | |
* @param h1 The first set of ijk coordinates. | |
* @param h2 The second set of ijk coordinates. | |
* @param sum The sum of the two sets of ijk coordinates. | |
*/ | |
void _ijkAdd(const CoordIJK *h1, const CoordIJK *h2, CoordIJK *sum) { | |
sum->i = h1->i + h2->i; | |
sum->j = h1->j + h2->j; | |
sum->k = h1->k + h2->k; | |
} | |
/** | |
* Subtract two ijk coordinates. | |
* | |
* @param h1 The first set of ijk coordinates. | |
* @param h2 The second set of ijk coordinates. | |
* @param diff The difference of the two sets of ijk coordinates (h1 - h2). | |
*/ | |
void _ijkSub(const CoordIJK *h1, const CoordIJK *h2, | |
CoordIJK *diff) { | |
diff->i = h1->i - h2->i; | |
diff->j = h1->j - h2->j; | |
diff->k = h1->k - h2->k; | |
} | |
/** | |
* Uniformly scale ijk coordinates by a scalar. Works in place. | |
* | |
* @param c The ijk coordinates to scale. | |
* @param factor The scaling factor. | |
*/ | |
void _ijkScale(CoordIJK *c, int factor) { | |
c->i *= factor; | |
c->j *= factor; | |
c->k *= factor; | |
} | |
/** | |
* Returns whether or not two ijk coordinates contain exactly the same | |
* component values. | |
* | |
* @param c1 The first set of ijk coordinates. | |
* @param c2 The second set of ijk coordinates. | |
* @return 1 if the two addresses match, 0 if they do not. | |
*/ | |
int _ijkMatches(const CoordIJK *c1, const CoordIJK *c2) { | |
return (c1->i == c2->i && c1->j == c2->j && c1->k == c2->k); | |
} | |
/** | |
* Find the normalized ijk coordinates of the indexing parent of a cell in a | |
* counter-clockwise aperture 7 grid. Works in place. | |
* | |
* @param ijk The ijk coordinates. | |
*/ | |
void _upAp7(CoordIJK *ijk) { | |
// convert to CoordIJ | |
int i = ijk->i - ijk->k; | |
int j = ijk->j - ijk->k; | |
ijk->i = (int)round((3 * i - j) / 7.0); | |
ijk->j = (int)round((i + 2 * j) / 7.0); | |
ijk->k = 0; | |
_ijkNormalize(ijk); | |
} | |
/** | |
* Find the normalized ijk coordinates of the indexing parent of a cell in a | |
* clockwise aperture 7 grid. Works in place. | |
* | |
* @param ijk The ijk coordinates. | |
*/ | |
void _upAp7r(CoordIJK *ijk) { | |
// convert to CoordIJ | |
int i = ijk->i - ijk->k; | |
int j = ijk->j - ijk->k; | |
ijk->i = (int)round((2 * i + j) / 7.0); | |
ijk->j = (int)round((3 * j - i) / 7.0); | |
ijk->k = 0; | |
_ijkNormalize(ijk); | |
} | |
/** | |
* Find the normalized ijk coordinates of the hex centered on the indicated | |
* hex at the next finer aperture 7 counter-clockwise resolution. Works in | |
* place. | |
* | |
* @param ijk The ijk coordinates. | |
*/ | |
void _downAp7(CoordIJK *ijk) { | |
// res r unit vectors in res r+1 | |
CoordIJK iVec = {3, 0, 1}; | |
CoordIJK jVec = {1, 3, 0}; | |
CoordIJK kVec = {0, 1, 3}; | |
_ijkScale(&iVec, ijk->i); | |
_ijkScale(&jVec, ijk->j); | |
_ijkScale(&kVec, ijk->k); | |
_ijkAdd(&iVec, &jVec, ijk); | |
_ijkAdd(ijk, &kVec, ijk); | |
_ijkNormalize(ijk); | |
} | |
/** | |
* Find the normalized ijk coordinates of the hex centered on the indicated | |
* hex at the next finer aperture 7 clockwise resolution. Works in place. | |
* | |
* @param ijk The ijk coordinates. | |
*/ | |
void _downAp7r(CoordIJK *ijk) { | |
// res r unit vectors in res r+1 | |
CoordIJK iVec = {3, 1, 0}; | |
CoordIJK jVec = {0, 3, 1}; | |
CoordIJK kVec = {1, 0, 3}; | |
_ijkScale(&iVec, ijk->i); | |
_ijkScale(&jVec, ijk->j); | |
_ijkScale(&kVec, ijk->k); | |
_ijkAdd(&iVec, &jVec, ijk); | |
_ijkAdd(ijk, &kVec, ijk); | |
_ijkNormalize(ijk); | |
} | |
// /** | |
// * Determines the H3 digit corresponding to a unit vector in ijk coordinates. | |
// * | |
// * @param ijk The ijk coordinates; must be a unit vector. | |
// * @return The H3 digit (0-6) corresponding to the ijk unit vector, or | |
// * INVALID_DIGIT on failure. | |
// */ | |
// Direction _unitIjkToDigit(const CoordIJK *ijk) { | |
// CoordIJK c = *ijk; | |
// _ijkNormalize(&c); | |
// Direction digit = INVALID_DIGIT; | |
// for (Direction i = CENTER_DIGIT; i < NUM_DIGITS; i++) { | |
// if (_ijkMatches(&c, &UNIT_VECS[i])) { | |
// digit = i; | |
// break; | |
// } | |
// } | |
// return digit; | |
// } | |
/** | |
* Determines the H3 digit corresponding to a unit vector in ijk coordinates. | |
* | |
* @param ijk The ijk coordinates; must be a unit vector. | |
* @return The H3 digit (0-6) corresponding to the ijk unit vector, or | |
* INVALID_DIGIT on failure. | |
*/ | |
Direction _unitIjkToDigit(const CoordIJK *ijk) { | |
CoordIJK c = *ijk; | |
_ijkNormalize(&c); | |
// TODO: Replicate enum loop | |
// CENTER_DIGIT = 0, | |
Direction digit = CENTER_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
// K_AXES_DIGIT = 1, | |
digit = K_AXES_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
// J_AXES_DIGIT = 2, | |
digit = J_AXES_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
// JK_AXES_DIGIT = J_AXES_DIGIT | K_AXES_DIGIT, /* 3 */ | |
digit = JK_AXES_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
// I_AXES_DIGIT = 4, | |
digit = I_AXES_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
// IK_AXES_DIGIT = I_AXES_DIGIT | K_AXES_DIGIT, /* 5 */ | |
digit = IK_AXES_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
// IJ_AXES_DIGIT = I_AXES_DIGIT | J_AXES_DIGIT, /* 6 */ | |
digit = IJ_AXES_DIGIT; | |
if (_ijkMatches(&c, &UNIT_VECS[digit])) { | |
return digit; | |
} | |
return INVALID_DIGIT; | |
} | |
/** @brief Find base cell given FaceIJK. | |
* | |
* Given the face number and a resolution 0 ijk+ coordinate in that face's | |
* face-centered ijk coordinate system, return the number of 60' ccw rotations | |
* to rotate into the coordinate system of the base cell at that coordinates. | |
* | |
* Valid ijk+ lookup coordinates are from (0, 0, 0) to (2, 2, 2). | |
*/ | |
int _faceIjkToBaseCellCCWrot60(const FaceIJK *h) { | |
return faceIjkBaseCells[h->face][h->coord.i][h->coord.j][h->coord.k].ccwRot60; | |
} | |
/** @brief Return whether or not the indicated base cell is a pentagon. */ | |
int _isBaseCellPentagon(int baseCell) { | |
if (baseCell < 0 || baseCell >= NUM_BASE_CELLS) { | |
// Base cells less than zero can not be represented in an index | |
return false; | |
} | |
return baseCellData[baseCell].isPentagon; | |
} | |
/** | |
* Returns the highest resolution non-zero digit in an H3Index. | |
* @param h The H3Index. | |
* @return The highest resolution non-zero digit in the H3Index. | |
*/ | |
Direction _h3LeadingNonZeroDigit(H3Index h) { | |
for (int r = 1; r <= H3_GET_RESOLUTION(h); r++) | |
if (H3_GET_INDEX_DIGIT(h, r)) | |
return H3_GET_INDEX_DIGIT(h, r); | |
// if we're here it's all 0's | |
return CENTER_DIGIT; | |
} | |
/** @brief Return whether or not the tested face is a cw offset face. | |
*/ | |
bool _baseCellIsCwOffset(int baseCell, int testFace) { | |
return baseCellData[baseCell].cwOffsetPent[0] == testFace || | |
baseCellData[baseCell].cwOffsetPent[1] == testFace; | |
} | |
/** | |
* Rotates indexing digit 60 degrees counter-clockwise. Returns result. | |
* | |
* @param digit Indexing digit (between 1 and 6 inclusive) | |
*/ | |
Direction _rotate60ccw(Direction digit) { | |
switch (digit) { | |
case K_AXES_DIGIT: | |
return IK_AXES_DIGIT; | |
case IK_AXES_DIGIT: | |
return I_AXES_DIGIT; | |
case I_AXES_DIGIT: | |
return IJ_AXES_DIGIT; | |
case IJ_AXES_DIGIT: | |
return J_AXES_DIGIT; | |
case J_AXES_DIGIT: | |
return JK_AXES_DIGIT; | |
case JK_AXES_DIGIT: | |
return K_AXES_DIGIT; | |
default: | |
return digit; | |
} | |
} | |
/** | |
* Rotates indexing digit 60 degrees clockwise. Returns result. | |
* | |
* @param digit Indexing digit (between 1 and 6 inclusive) | |
*/ | |
Direction _rotate60cw(Direction digit) { | |
switch (digit) { | |
case K_AXES_DIGIT: | |
return JK_AXES_DIGIT; | |
case JK_AXES_DIGIT: | |
return J_AXES_DIGIT; | |
case J_AXES_DIGIT: | |
return IJ_AXES_DIGIT; | |
case IJ_AXES_DIGIT: | |
return I_AXES_DIGIT; | |
case I_AXES_DIGIT: | |
return IK_AXES_DIGIT; | |
case IK_AXES_DIGIT: | |
return K_AXES_DIGIT; | |
default: | |
return digit; | |
} | |
} | |
/** | |
* Rotate an H3Index 60 degrees counter-clockwise. | |
* @param h The H3Index. | |
*/ | |
H3Index _h3Rotate60ccw(H3Index h) { | |
for (int r = 1, res = H3_GET_RESOLUTION(h); r <= res; r++) { | |
Direction oldDigit = H3_GET_INDEX_DIGIT(h, r); | |
H3_SET_INDEX_DIGIT(h, r, _rotate60ccw(oldDigit)); | |
} | |
return h; | |
} | |
/** | |
* Rotate an H3Index 60 degrees clockwise. | |
* @param h The H3Index. | |
*/ | |
H3Index _h3Rotate60cw(H3Index h) { | |
for (int r = 1, res = H3_GET_RESOLUTION(h); r <= res; r++) { | |
H3_SET_INDEX_DIGIT(h, r, _rotate60cw(H3_GET_INDEX_DIGIT(h, r))); | |
} | |
return h; | |
} | |
/** | |
* Rotate an H3Index 60 degrees counter-clockwise about a pentagonal center. | |
* @param h The H3Index. | |
*/ | |
H3Index _h3RotatePent60ccw(H3Index h) { | |
// rotate in place; skips any leading 1 digits (k-axis) | |
int foundFirstNonZeroDigit = 0; | |
for (int r = 1, res = H3_GET_RESOLUTION(h); r <= res; r++) { | |
// rotate this digit | |
H3_SET_INDEX_DIGIT(h, r, _rotate60ccw(H3_GET_INDEX_DIGIT(h, r))); | |
// look for the first non-zero digit so we | |
// can adjust for deleted k-axes sequence | |
// if necessary | |
if (!foundFirstNonZeroDigit && H3_GET_INDEX_DIGIT(h, r) != 0) { | |
foundFirstNonZeroDigit = 1; | |
// adjust for deleted k-axes sequence | |
if (_h3LeadingNonZeroDigit(h) == K_AXES_DIGIT) | |
h = _h3Rotate60ccw(h); | |
} | |
} | |
return h; | |
} | |
/** | |
* Convert an FaceIJK address to the corresponding H3Index. | |
* @param fijk The FaceIJK address. | |
* @param res The cell resolution. | |
* @return The encoded H3Index (or H3_NULL on failure). | |
*/ | |
H3Index _faceIjkToH3(const FaceIJK *fijk, int res) { | |
// initialize the index | |
H3Index h = H3_INIT; | |
H3_SET_MODE(h, H3_CELL_MODE); | |
H3_SET_RESOLUTION(h, res); | |
// check for res 0/base cell | |
if (res == 0) { | |
if (fijk->coord.i > MAX_FACE_COORD || fijk->coord.j > MAX_FACE_COORD || | |
fijk->coord.k > MAX_FACE_COORD) { | |
// out of range input | |
return H3_NULL; | |
} | |
H3_SET_BASE_CELL(h, _faceIjkToBaseCell(fijk)); | |
return h; | |
} | |
// we need to find the correct base cell FaceIJK for this H3 index; | |
// start with the passed in face and resolution res ijk coordinates | |
// in that face's coordinate system | |
FaceIJK fijkBC = *fijk; | |
// return h; | |
// build the H3Index from finest res up | |
// adjust r for the fact that the res 0 base cell offsets the indexing | |
// digits | |
CoordIJK *ijk = &fijkBC.coord; | |
for (int r = res - 1; r >= 0; r--) { | |
CoordIJK lastIJK = *ijk; | |
CoordIJK lastCenter; | |
if (isResolutionClassIII(r + 1)) { | |
// rotate ccw | |
_upAp7(ijk); | |
lastCenter = *ijk; | |
_downAp7(&lastCenter); | |
} else { | |
// rotate cw | |
_upAp7r(ijk); | |
lastCenter = *ijk; | |
_downAp7r(&lastCenter); | |
} | |
CoordIJK diff; | |
_ijkSub(&lastIJK, &lastCenter, &diff); | |
_ijkNormalize(&diff); | |
H3_SET_INDEX_DIGIT(h, r + 1, _unitIjkToDigit(&diff)); | |
} | |
// fijkBC should now hold the IJK of the base cell in the | |
// coordinate system of the current face | |
if (fijkBC.coord.i > MAX_FACE_COORD || fijkBC.coord.j > MAX_FACE_COORD || | |
fijkBC.coord.k > MAX_FACE_COORD) { | |
// out of range input | |
return H3_NULL; | |
} | |
// lookup the correct base cell | |
int baseCell = _faceIjkToBaseCell(&fijkBC); | |
H3_SET_BASE_CELL(h, baseCell); | |
// rotate if necessary to get canonical base cell orientation | |
// for this base cell | |
int numRots = _faceIjkToBaseCellCCWrot60(&fijkBC); | |
if (_isBaseCellPentagon(baseCell)) { | |
// force rotation out of missing k-axes sub-sequence | |
if (_h3LeadingNonZeroDigit(h) == K_AXES_DIGIT) { | |
// check for a cw/ccw offset face; default is ccw | |
if (_baseCellIsCwOffset(baseCell, fijkBC.face)) { | |
h = _h3Rotate60cw(h); | |
} else { | |
h = _h3Rotate60ccw(h); | |
} | |
} | |
for (int i = 0; i < numRots; i++) | |
h = _h3RotatePent60ccw(h); | |
} else { | |
for (int i = 0; i < numRots; i++) { | |
h = _h3Rotate60ccw(h); | |
} | |
} | |
return h; | |
} | |
/** | |
* Encodes a coordinate on the sphere to the H3 index of the containing cell at | |
* the specified resolution. | |
* | |
* Returns 0 on invalid input. | |
* | |
* @param g The spherical coordinates to encode. | |
* @param res The desired H3 resolution for the encoding. | |
* @param out The encoded H3Index. | |
* @returns E_SUCCESS (0) on success, another value otherwise | |
*/ | |
H3Error latLngToCell(const LatLng *g, int res, H3Index *out) { | |
if (res < 0 || res > MAX_H3_RES) { | |
return E_RES_DOMAIN; | |
} | |
if (!isfinite(g->lat) || !isfinite(g->lng)) { | |
return E_LATLNG_DOMAIN; | |
} | |
FaceIJK fijk; | |
_geoToFaceIjk(g, res, &fijk); | |
*out = _faceIjkToH3(&fijk, res); | |
if (*out) { | |
return E_SUCCESS; | |
} else { | |
return E_FAILED; | |
} | |
} | |
__kernel void h3_paired( | |
__global const float *lat_g, __global const float *lng_g, const unsigned int resolution, __global ulong *res_g) | |
{ | |
int gid = get_global_id(0); | |
LatLng location; | |
H3Index indexed; | |
location.lat = degsToRads(lat_g[gid]); | |
location.lng = degsToRads(lng_g[gid]); | |
if (latLngToCell(&location, resolution, &indexed) != E_SUCCESS) { | |
res_g[gid] = H3_NULL; | |
} else { | |
res_g[gid] = indexed; | |
} | |
} | |
__kernel void h3_single( | |
__global const float *pts_g, const unsigned int resolution, __global ulong *res_g) | |
{ | |
int gid = get_global_id(0); | |
int idx = gid * 2; | |
LatLng location; | |
H3Index indexed; | |
location.lat = degsToRads(pts_g[idx]); | |
location.lng = degsToRads(pts_g[idx+1]); | |
if (latLngToCell(&location, resolution, &indexed) != E_SUCCESS) { | |
res_g[gid] = H3_NULL; | |
} else { | |
res_g[gid] = indexed; | |
} | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# python 3.10.4 on pop_os 20.04 (debian). pyopencl==2022.1.3 numpy=1.22.3 h3==3.7.3 | |
from time import time | |
import warnings | |
import numpy as np | |
import pyopencl as cl | |
from functools import lru_cache | |
with warnings.catch_warnings(): | |
warnings.filterwarnings("ignore", category=UserWarning) | |
import h3.unstable.vect as h3_vect # h3 helpfully and correctly notes this is unstable | |
def generate_lat_lon_points(count: int = 50, seed: int = 0): | |
rng = np.random.default_rng(seed) | |
lat = ((rng.random(count) * 2 - 1) * 90.0).astype(np.float32) | |
lng = ((rng.random(count) * 2 - 1) * 180.0).astype(np.float32) | |
return lat, lng | |
@lru_cache | |
def get_ctx(): | |
try: | |
platforms = cl.get_platforms() | |
except Exception as why: | |
# logger.error("Unable to get OpenCL platform") | |
raise why | |
# if len(platforms) > 1: | |
# logger.warning(f"Detected {len(platforms)} OpenCL platforms. Defaulting to {platforms[0]}") | |
ctx = cl.Context(dev_type=cl.device_type.ALL,properties=[(cl.context_properties.PLATFORM, platforms[0])]) | |
return ctx | |
@lru_cache | |
def get_cl_h3(): | |
with open("h3.cl", 'r') as fp: | |
h3_code = fp.read() | |
return h3_code | |
def opencl(lat_np:np.ndarray, lng_np:np.ndarray, resolution:int)->np.ndarray: | |
assert lat_np.dtype == np.float32 | |
assert lng_np.dtype == np.float32 | |
ctx = get_ctx() | |
queue = cl.CommandQueue(ctx) | |
prg = cl.Program(ctx, get_cl_h3()).build() | |
mf = cl.mem_flags | |
lat_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=lat_np) | |
lng_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=lng_np) | |
res_np = np.zeros((lat_np.shape[0]), dtype=np.uint64) | |
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, res_np.nbytes) | |
h3_paired_knl = prg.h3_paired # Use this Kernel object for repeated calls | |
h3_paired_knl(queue, (lat_np.shape), None, lat_g, lng_g, np.int32(resolution), res_g) | |
cl.enqueue_copy(queue, res_np, res_g) | |
return res_np | |
def trail(lat_np, lng_np, testing): | |
resolution = 4 | |
start = time() | |
res_np = opencl(lat_np, lng_np, resolution) | |
print(f"OpenCL for {lat_np.shape[0]:,} in {time() - start: >6.5f} seconds") | |
if testing: | |
expected = h3_vect.geo_to_h3(lat_np, lng_np, resolution) | |
np.testing.assert_allclose(expected, res_np) | |
def trail_alpha(lat_np, lng_np, testing): | |
trail(lat_np, lng_np, testing) | |
def trail_beta(lat_np, lng_np, testing): | |
trail(lat_np[:-1], lng_np[:-1], testing) | |
def main(): | |
count = 18_001_933 | |
testing = False # enable to verify results | |
lat_np, lng_np = generate_lat_lon_points(count, seed=8) | |
trail_alpha(lat_np, lng_np, testing) | |
trail_beta(lat_np, lng_np, testing) | |
trail_alpha(lat_np, lng_np, testing) | |
if __name__ == "__main__": | |
main() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment