Skip to content

Instantly share code, notes, and snippets.

@user01

user01/clinfo.txt Secret

Created May 9, 2022 02:42
Show Gist options
  • Save user01/514a9ad748858310783c297d42659c45 to your computer and use it in GitHub Desktop.
Save user01/514a9ad748858310783c297d42659c45 to your computer and use it in GitHub Desktop.
Download ZIP
pyopencl h3 example
Raw
clinfo.txt
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
Raw
h3.cl
/** 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;
}
}
Raw
pycl_exp.py
# 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