NanoVDB.h

Go to the documentation of this file.

// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0

/*!
    \file   NanoVDB.h

    \author Ken Museth

    \date  January 8, 2020

    \brief Implements a light-weight self-contained VDB data-structure in a
           single file! In other words, this is a significantly watered-down
           version of the OpenVDB implementation, with few dependencies - so
           a one-stop-shop for a minimalistic VDB data structure that run on
           most platforms!

    \note It is important to note that NanoVDB (by design) is a read-only
          sparse GPU (and CPU) friendly data structure intended for applications
          like rendering and collision detection. As such it obviously lacks
          a lot of the functionality and features of OpenVDB grids. NanoVDB
          is essentially a compact linearized (or serialized) representation of
          an OpenVDB tree with getValue methods only. For best performance use
          the ReadAccessor::getValue method as opposed to the Tree::getValue
          method. Note that since a ReadAccessor caches previous access patterns
          it is by design not thread-safe, so use one instantiation per thread
          (it is very light-weight). Also, it is not safe to copy accessors between
          the GPU and CPU! In fact, client code should only interface
          with the API of the Grid class (all other nodes of the NanoVDB data
          structure can safely be ignored by most client codes)!


    \warning NanoVDB grids can only be constructed via tools like createNanoGrid
             or the GridBuilder. This explains why none of the grid nodes defined below
             have public constructors or destructors.

    \details Please see the following paper for more details on the data structure:
          K. Museth, “VDB: High-Resolution Sparse Volumes with Dynamic Topology”,
          ACM Transactions on Graphics 32(3), 2013, which can be found here:
          http://www.museth.org/Ken/Publications_files/Museth_TOG13.pdf

          NanoVDB was first published there: https://dl.acm.org/doi/fullHtml/10.1145/3450623.3464653


    Overview: This file implements the following fundamental class that when combined
          forms the backbone of the VDB tree data structure:

          Coord- a signed integer coordinate
          Vec3 - a 3D vector
          Vec4 - a 4D vector
          BBox - a bounding box
          Mask - a bitmask essential to the non-root tree nodes
          Map  - an affine coordinate transformation
          Grid - contains a Tree and a map for world<->index transformations. Use
                 this class as the main API with client code!
          Tree - contains a RootNode and getValue methods that should only be used for debugging
          RootNode - the top-level node of the VDB data structure
          InternalNode - the internal nodes of the VDB data structure
          LeafNode - the lowest level tree nodes that encode voxel values and state
          ReadAccessor - implements accelerated random access operations

    Semantics: A VDB data structure encodes values and (binary) states associated with
          signed integer coordinates. Values encoded at the leaf node level are
          denoted voxel values, and values associated with other tree nodes are referred
          to as tile values, which by design cover a larger coordinate index domain.


    Memory layout:

    It's important to emphasize that all the grid data (defined below) are explicitly 32 byte
    aligned, which implies that any memory buffer that contains a NanoVDB grid must also be at
    32 byte aligned. That is, the memory address of the beginning of a buffer (see ascii diagram below)
    must be divisible by 32, i.e. uintptr_t(&buffer)%32 == 0! If this is not the case, the C++ standard
    says the behaviour is undefined! Normally this is not a concerns on GPUs, because they use 256 byte
    aligned allocations, but the same cannot be said about the CPU.

    GridData is always at the very beginning of the buffer immediately followed by TreeData!
    The remaining nodes and blind-data are allowed to be scattered throughout the buffer,
    though in practice they are arranged as:

    GridData: 672 bytes (e.g. magic, checksum, major, flags, index, count, size, name, map, world bbox, voxel size, class, type, offset, count)

    TreeData: 64 bytes (node counts and byte offsets)

    ... optional padding ...

    RootData: size depends on ValueType (index bbox, voxel count, tile count, min/max/avg/standard deviation)

    Array of: RootData::Tile

    ... optional padding ...

    Array of: Upper InternalNodes of size 32^3:  bbox, two bit masks, 32768 tile values, and min/max/avg/standard deviation values

    ... optional padding ...

    Array of: Lower InternalNodes of size 16^3:  bbox, two bit masks, 4096 tile values, and min/max/avg/standard deviation values

    ... optional padding ...

    Array of: LeafNodes of size 8^3: bbox, bit masks, 512 voxel values, and min/max/avg/standard deviation values


    Notation: "]---[" implies it has optional padding, and "][" implies zero padding

    [GridData(672B)][TreeData(64B)]---[RootData][N x Root::Tile]---[InternalData<5>]---[InternalData<4>]---[LeafData<3>]---[BLINDMETA...]---[BLIND0]---[BLIND1]---etc.
    ^                                 ^         ^                  ^                   ^                   ^
    |                                 |         |                  |                   |                   |
    +-- Start of 32B aligned buffer   |         |                  |                   |                   +-- Node0::DataType* leafData
        GridType::DataType* gridData  |         |                  |                   |
                                      |         |                  |                   +-- Node1::DataType* lowerData
       RootType::DataType* rootData --+         |                  |
                                                |                  +-- Node2::DataType* upperData
                                                |
                                                +-- RootType::DataType::Tile* tile

*/

#ifndef NANOVDB_NANOVDB_H_HAS_BEEN_INCLUDED
#define NANOVDB_NANOVDB_H_HAS_BEEN_INCLUDED

// NANOVDB_MAGIC_NUMBER is currently used for both grids and files (starting with v32.6.0)
// NANOVDB_MAGIC_GRID will soon be used exclusively for grids
// NANOVDB_MAGIC_FILE will soon be used exclusively for files
// NANOVDB_MAGIC_NODE will soon be used exclusively for NodeManager
//                              | : 0 in 30 corresponds to 0 in NanoVDB0
#define NANOVDB_MAGIC_NUMBER 0x304244566f6e614eUL // "NanoVDB0" in hex - little endian (uint64_t)
#define NANOVDB_MAGIC_GRID   0x314244566f6e614eUL // "NanoVDB1" in hex - little endian (uint64_t)
#define NANOVDB_MAGIC_FILE   0x324244566f6e614eUL // "NanoVDB2" in hex - little endian (uint64_t)
#define NANOVDB_MAGIC_NODE   0x334244566f6e614eUL // "NanoVDB3" in hex - little endian (uint64_t)
#define NANOVDB_MAGIC_MASK   0x00FFFFFFFFFFFFFFUL // use this mask to remove the number
//#define NANOVDB_USE_NEW_MAGIC_NUMBERS// used to enable use of the new magic numbers described above

#define NANOVDB_MAJOR_VERSION_NUMBER 32 // reflects changes to the ABI and hence also the file format
#define NANOVDB_MINOR_VERSION_NUMBER 6 //  reflects changes to the API but not ABI
#define NANOVDB_PATCH_VERSION_NUMBER 0 //  reflects changes that does not affect the ABI or API

#define TBB_SUPPRESS_DEPRECATED_MESSAGES 1

// This replaces a Coord key at the root level with a single uint64_t
#define NANOVDB_USE_SINGLE_ROOT_KEY

// This replaces three levels of Coord keys in the ReadAccessor with one Coord
//#define NANOVDB_USE_SINGLE_ACCESSOR_KEY

// Use this to switch between std::ofstream or FILE implementations
//#define NANOVDB_USE_IOSTREAMS

// Use this to switch between old and new accessor methods
#define NANOVDB_NEW_ACCESSOR_METHODS

#define NANOVDB_FPN_BRANCHLESS

// Do not change this value! 32 byte alignment is fixed in NanoVDB
#define NANOVDB_DATA_ALIGNMENT 32

#if !defined(NANOVDB_ALIGN)
#define NANOVDB_ALIGN(n) alignas(n)
#endif // !defined(NANOVDB_ALIGN)

#ifdef __CUDACC_RTC__

typedef signed char        int8_t;
typedef short              int16_t;
typedef int                int32_t;
typedef long long          int64_t;
typedef unsigned char      uint8_t;
typedef unsigned int       uint32_t;
typedef unsigned short     uint16_t;
typedef unsigned long long uint64_t;

#define NANOVDB_ASSERT(x)

#define UINT64_C(x)  (x ## ULL)

#else // !__CUDACC_RTC__

#include <stdlib.h> //    for abs in clang7
#include <stdint.h> //    for types like int32_t etc
#include <stddef.h> //    for size_t type
#include <cassert> //     for assert
#include <cstdio> //      for snprintf
#include <cmath> //       for sqrt and fma
#include <limits> //      for numeric_limits
#include <utility>//      for std::move
#ifdef NANOVDB_USE_IOSTREAMS
#include <fstream>//      for read/writeUncompressedGrids
#endif
// All asserts can be disabled here, even for debug builds
#if 1
#define NANOVDB_ASSERT(x) assert(x)
#else
#define NANOVDB_ASSERT(x)
#endif

#if defined(NANOVDB_USE_INTRINSICS) && defined(_MSC_VER)
#include <intrin.h>
#pragma intrinsic(_BitScanReverse)
#pragma intrinsic(_BitScanForward)
#pragma intrinsic(_BitScanReverse64)
#pragma intrinsic(_BitScanForward64)
#endif

#endif // __CUDACC_RTC__

#if defined(__CUDACC__) || defined(__HIP__)
// Only define __hostdev__ when using NVIDIA CUDA or HIP compilers
#ifndef __hostdev__
#define __hostdev__ __host__ __device__ // Runs on the CPU and GPU, called from the CPU or the GPU
#endif
#else
// Dummy definitions of macros only defined by CUDA and HIP compilers
#ifndef __hostdev__
#define __hostdev__ // Runs on the CPU and GPU, called from the CPU or the GPU
#endif
#ifndef __global__
#define __global__ // Runs on the GPU, called from the CPU or the GPU
#endif
#ifndef __device__
#define __device__ // Runs on the GPU, called from the GPU
#endif
#ifndef __host__
#define __host__ // Runs on the CPU, called from the CPU
#endif

#endif // if defined(__CUDACC__) || defined(__HIP__)

// The following macro will suppress annoying warnings when nvcc
// compiles functions that call (host) intrinsics (which is perfectly valid)
#if defined(_MSC_VER) && defined(__CUDACC__)
#define NANOVDB_HOSTDEV_DISABLE_WARNING __pragma("hd_warning_disable")
#elif defined(__GNUC__) && defined(__CUDACC__)
#define NANOVDB_HOSTDEV_DISABLE_WARNING _Pragma("hd_warning_disable")
#else
#define NANOVDB_HOSTDEV_DISABLE_WARNING
#endif

// Define compiler warnings that work with all compilers
//#if defined(_MSC_VER)
//#define NANO_WARNING(msg) _pragma("message" #msg)
//#else
//#define NANO_WARNING(msg) _Pragma("message" #msg)
//#endif

// A portable implementation of offsetof - unfortunately it doesn't work with static_assert
#define NANOVDB_OFFSETOF(CLASS, MEMBER) ((int)(size_t)((char*)&((CLASS*)0)->MEMBER - (char*)0))

namespace nanovdb {

// --------------------------> Build types <------------------------------------

/// @brief Dummy type for a voxel whose value equals an offset into an external value array
class ValueIndex{};

/// @brief Dummy type for a voxel whose value equals an offset into an external value array of active values
class ValueOnIndex{};

/// @brief Like @c ValueIndex but with a mutable mask
class ValueIndexMask{};

/// @brief Like @c ValueOnIndex but with a mutable mask
class ValueOnIndexMask{};

/// @brief Dummy type for a voxel whose value equals its binary active state
class ValueMask{};

/// @brief Dummy type for a 16 bit floating point values (placeholder for IEEE 754 Half)
class Half{};

/// @brief Dummy type for a 4bit quantization of float point values
class Fp4{};

/// @brief Dummy type for a 8bit quantization of float point values
class Fp8{};

/// @brief Dummy type for a 16bit quantization of float point values
class Fp16{};

/// @brief Dummy type for a variable bit quantization of floating point values
class FpN{};

/// @brief Dummy type for indexing points into voxels
class Point{};

// --------------------------> GridType <------------------------------------

/// @brief List of types that are currently supported by NanoVDB
///
/// @note To expand on this list do:
///       1) Add the new type between Unknown and End in the enum below
///       2) Add the new type to OpenToNanoVDB::processGrid that maps OpenVDB types to GridType
///       3) Verify that the ConvertTrait in NanoToOpenVDB.h works correctly with the new type
///       4) Add the new type to mapToGridType (defined below) that maps NanoVDB types to GridType
///       5) Add the new type to toStr (defined below)
enum class GridType : uint32_t { Unknown = 0, //  unknown value type - should rarely be used
                                 Float = 1, //  single precision floating point value
                                 Double = 2, //  double precision floating point value
                                 Int16 = 3, //  half precision signed integer value
                                 Int32 = 4, //  single precision signed integer value
                                 Int64 = 5, //  double precision signed integer value
                                 Vec3f = 6, //  single precision floating 3D vector
                                 Vec3d = 7, //  double precision floating 3D vector
                                 Mask = 8, //  no value, just the active state
                                 Half = 9, //  half precision floating point value (placeholder for IEEE 754 Half)
                                 UInt32 = 10, // single precision unsigned integer value
                                 Boolean = 11, // boolean value, encoded in bit array
                                 RGBA8 = 12, // RGBA packed into 32bit word in reverse-order, i.e. R is lowest byte.
                                 Fp4 = 13, // 4bit quantization of floating point value
                                 Fp8 = 14, // 8bit quantization of floating point value
                                 Fp16 = 15, // 16bit quantization of floating point value
                                 FpN = 16, // variable bit quantization of floating point value
                                 Vec4f = 17, // single precision floating 4D vector
                                 Vec4d = 18, // double precision floating 4D vector
                                 Index = 19, // index into an external array of active and inactive values
                                 OnIndex = 20, // index into an external array of active values
                                 IndexMask = 21, // like Index but with a mutable mask
                                 OnIndexMask = 22, // like OnIndex but with a mutable mask
                                 PointIndex = 23, // voxels encode indices to co-located points
                                 Vec3u8 = 24, // 8bit quantization of floating point 3D vector (only as blind data)
                                 Vec3u16 = 25, // 16bit quantization of floating point 3D vector (only as blind data)
                                 End = 26 }; // should never be used

#ifndef __CUDACC_RTC__
/// @brief Maps a GridType to a c-string
/// @param gridType GridType to be mapped to a string
/// @return Retuns a c-string used to describe a GridType
inline const char* toStr(GridType gridType)
{
    static const char* LUT[] = {"?", "float", "double", "int16", "int32", "int64", "Vec3f", "Vec3d", "Mask", "Half",
                                "uint32", "bool", "RGBA8", "Float4", "Float8", "Float16", "FloatN", "Vec4f", "Vec4d",
                                "Index", "OnIndex", "IndexMask", "OnIndexMask", "PointIndex", "Vec3u8", "Vec3u16", "End"};
    static_assert(sizeof(LUT) / sizeof(char*) - 1 == int(GridType::End), "Unexpected size of LUT");
    return LUT[static_cast<int>(gridType)];
}
#endif

// --------------------------> GridClass <------------------------------------

/// @brief Classes (superset of OpenVDB) that are currently supported by NanoVDB
enum class GridClass : uint32_t { Unknown = 0,
                                  LevelSet = 1, // narrow band level set, e.g. SDF
                                  FogVolume = 2, // fog volume, e.g. density
                                  Staggered = 3, // staggered MAC grid, e.g. velocity
                                  PointIndex = 4, // point index grid
                                  PointData = 5, // point data grid
                                  Topology = 6, // grid with active states only (no values)
                                  VoxelVolume = 7, // volume of geometric cubes, e.g. colors cubes in Minecraft
                                  IndexGrid = 8, // grid whose values are offsets, e.g. into an external array
                                  TensorGrid = 9, // Index grid for indexing learnable tensor features
                                  End = 10 };

#ifndef __CUDACC_RTC__
/// @brief Retuns a c-string used to describe a GridClass
inline const char* toStr(GridClass gridClass)
{
    static const char* LUT[] = {"?", "SDF", "FOG", "MAC", "PNTIDX", "PNTDAT", "TOPO", "VOX", "INDEX", "TENSOR", "END"};
    static_assert(sizeof(LUT) / sizeof(char*) - 1 == int(GridClass::End), "Unexpected size of LUT");
    return LUT[static_cast<int>(gridClass)];
}
#endif

// --------------------------> GridFlags <------------------------------------

/// @brief Grid flags which indicate what extra information is present in the grid buffer.
enum class GridFlags : uint32_t {
    HasLongGridName = 1 << 0, // grid name is longer than 256 characters
    HasBBox = 1 << 1, // nodes contain bounding-boxes of active values
    HasMinMax = 1 << 2, // nodes contain min/max of active values
    HasAverage = 1 << 3, // nodes contain averages of active values
    HasStdDeviation = 1 << 4, // nodes contain standard deviations of active values
    IsBreadthFirst = 1 << 5, // nodes are typically arranged breadth-first in memory
    End = 1 << 6, // use End - 1 as a mask for the 5 lower bit flags
};

#ifndef __CUDACC_RTC__
/// @brief Retuns a c-string used to describe a GridFlags
inline const char* toStr(GridFlags gridFlags)
{
    static const char* LUT[] = {"has long grid name",
                                "has bbox",
                                "has min/max",
                                "has average",
                                "has standard deviation",
                                "is breadth-first",
                                "end"};
    static_assert(1 << (sizeof(LUT) / sizeof(char*) - 1) == int(GridFlags::End), "Unexpected size of LUT");
    return LUT[static_cast<int>(gridFlags)];
}
#endif

// --------------------------> GridBlindData enums <------------------------------------

/// @brief Blind-data Classes that are currently supported by NanoVDB
enum class GridBlindDataClass : uint32_t { Unknown = 0,
                                           IndexArray = 1,
                                           AttributeArray = 2,
                                           GridName = 3,
                                           ChannelArray = 4,
                                           End = 5 };

/// @brief Blind-data Semantics that are currently understood by NanoVDB
enum class GridBlindDataSemantic : uint32_t { Unknown = 0,
                                              PointPosition = 1, // 3D coordinates in an unknown space
                                              PointColor = 2,
                                              PointNormal = 3,
                                              PointRadius = 4,
                                              PointVelocity = 5,
                                              PointId = 6,
                                              WorldCoords = 7, // 3D coordinates in world space, e.g. (0.056, 0.8, 1,8)
                                              GridCoords = 8, // 3D coordinates in grid space, e.g. (1.2, 4.0, 5.7), aka index-space
                                              VoxelCoords = 9, // 3D coordinates in voxel space, e.g. (0.2, 0.0, 0.7)
                                              End = 10 };

// --------------------------> is_same <------------------------------------

/// @brief C++11 implementation of std::is_same
/// @note When more than two arguments are provided value = T0==T1 || T0==T2 || ...
template<typename T0, typename T1, typename ...T>
struct is_same
{
    static constexpr bool value = is_same<T0, T1>::value || is_same<T0, T...>::value;
};

template<typename T0, typename T1>
struct is_same<T0, T1>
{
    static constexpr bool value = false;
};

template<typename T>
struct is_same<T, T>
{
    static constexpr bool value = true;
};

// --------------------------> is_floating_point <------------------------------------

/// @brief C++11 implementation of std::is_floating_point
template<typename T>
struct is_floating_point
{
    static constexpr bool value = is_same<T, float, double>::value;
};

// --------------------------> BuildTraits <------------------------------------

/// @brief Define static boolean tests for template build types
template<typename T>
struct BuildTraits
{
    // check if T is an index type
    static constexpr bool is_index     = is_same<T, ValueIndex, ValueIndexMask, ValueOnIndex, ValueOnIndexMask>::value;
    static constexpr bool is_onindex   = is_same<T, ValueOnIndex, ValueOnIndexMask>::value;
    static constexpr bool is_offindex  = is_same<T, ValueIndex, ValueIndexMask>::value;
    static constexpr bool is_indexmask = is_same<T, ValueIndexMask, ValueOnIndexMask>::value;
    // check if T is a compressed float type with fixed bit precision
    static constexpr bool is_FpX = is_same<T, Fp4, Fp8, Fp16>::value;
    // check if T is a compressed float type with fixed or variable bit precision
    static constexpr bool is_Fp = is_same<T, Fp4, Fp8, Fp16, FpN>::value;
    // check if T is a POD float type, i.e float or double
    static constexpr bool is_float = is_floating_point<T>::value;
    // check if T is a template specialization of LeafData<T>, i.e. has T mValues[512]
    static constexpr bool is_special = is_index || is_Fp || is_same<T, Point, bool, ValueMask>::value;
}; // BuildTraits

// --------------------------> enable_if <------------------------------------

/// @brief C++11 implementation of std::enable_if
template <bool, typename T = void>
struct enable_if
{
};

template <typename T>
struct enable_if<true, T>
{
    using type = T;
};

// --------------------------> disable_if <------------------------------------

template<bool, typename T = void>
struct disable_if
{
    typedef T type;
};

template<typename T>
struct disable_if<true, T>
{
};

// --------------------------> is_const <------------------------------------

template<typename T>
struct is_const
{
    static constexpr bool value = false;
};

template<typename T>
struct is_const<const T>
{
    static constexpr bool value = true;
};

// --------------------------> is_pointer <------------------------------------

/// @brief Trait used to identify template parameter that are pointers
/// @tparam T Template parameter to be tested
template<class T>
struct is_pointer
{
    static constexpr bool value = false;
};

/// @brief Template specialization of non-const pointers
/// @tparam T Template parameter to be tested
template<class T>
struct is_pointer<T*>
{
    static constexpr bool value = true;
};

/// @brief Template specialization of const pointers
/// @tparam T Template parameter to be tested
template<class T>
struct is_pointer<const T*>
{
    static constexpr bool value = true;
};

// --------------------------> remove_const <------------------------------------

/// @brief Trait use to const from type. Default implementation is just a pass-through
/// @tparam T Type
/// @details remove_pointer<float>::type = float
template<typename T>
struct remove_const
{
    using type = T;
};

/// @brief Template specialization of trait class use to remove const qualifier type from a type
/// @tparam T Type of the const type
/// @details remove_pointer<const float>::type = float
template<typename T>
struct remove_const<const T>
{
    using type = T;
};

// --------------------------> remove_reference <------------------------------------

/// @brief Trait use to remove reference, i.e. "&", qualifier from a type. Default implementation is just a pass-through
/// @tparam T Type
/// @details remove_pointer<float>::type = float
template <typename T>
struct remove_reference {using type = T;};

/// @brief Template specialization of trait class use to remove reference, i.e. "&", qualifier from a type
/// @tparam T Type of the reference
/// @details remove_pointer<float&>::type = float
template <typename T>
struct remove_reference<T&> {using type = T;};

// --------------------------> remove_pointer <------------------------------------

/// @brief Trait use to remove pointer, i.e. "*", qualifier from a type. Default implementation is just a pass-through
/// @tparam T Type
/// @details remove_pointer<float>::type = float
template <typename T>
struct remove_pointer {using type = T;};

/// @brief Template specialization of trait class use to to remove pointer, i.e. "*", qualifier from a type
/// @tparam T Type of the pointer
/// @details remove_pointer<float*>::type = float
template <typename T>
struct remove_pointer<T*> {using type = T;};

// --------------------------> match_const <------------------------------------

/// @brief Trait used to transfer the const-ness of a reference type to another type
/// @tparam T Type whose const-ness needs to match the reference type
/// @tparam ReferenceT Reference type that is not const
/// @details match_const<const int, float>::type = int
///          match_const<int, float>::type = int
template<typename T, typename ReferenceT>
struct match_const
{
    using type = typename remove_const<T>::type;
};

/// @brief Template specialization used to transfer the const-ness of a reference type to another type
/// @tparam T Type that will adopt the const-ness of the reference type
/// @tparam ReferenceT Reference type that is const
/// @details match_const<const int, const float>::type = const int
///          match_const<int, const float>::type = const int
template<typename T, typename ReferenceT>
struct match_const<T, const ReferenceT>
{
    using type = const typename remove_const<T>::type;
};

// --------------------------> is_specialization <------------------------------------

/// @brief Metafunction used to determine if the first template
///        parameter is a specialization of the class template
///        given in the second template parameter.
///
/// @details is_specialization<Vec3<float>, Vec3>::value == true;
///          is_specialization<Vec3f, Vec3>::value == true;
///          is_specialization<std::vector<float>, std::vector>::value == true;
template<typename AnyType, template<typename...> class TemplateType>
struct is_specialization
{
    static const bool value = false;
};
template<typename... Args, template<typename...> class TemplateType>
struct is_specialization<TemplateType<Args...>, TemplateType>
{
    static const bool value = true;
};

// --------------------------> BuildToValueMap <------------------------------------

/// @brief Maps one type (e.g. the build types above) to other (actual) types
template<typename T>
struct BuildToValueMap
{
    using Type = T;
    using type = T;
};

template<>
struct BuildToValueMap<ValueIndex>
{
    using Type = uint64_t;
    using type = uint64_t;
};

template<>
struct BuildToValueMap<ValueOnIndex>
{
    using Type = uint64_t;
    using type = uint64_t;
};

template<>
struct BuildToValueMap<ValueIndexMask>
{
    using Type = uint64_t;
    using type = uint64_t;
};

template<>
struct BuildToValueMap<ValueOnIndexMask>
{
    using Type = uint64_t;
    using type = uint64_t;
};

template<>
struct BuildToValueMap<ValueMask>
{
    using Type = bool;
    using type = bool;
};

template<>
struct BuildToValueMap<Half>
{
    using Type = float;
    using type = float;
};

template<>
struct BuildToValueMap<Fp4>
{
    using Type = float;
    using type = float;
};

template<>
struct BuildToValueMap<Fp8>
{
    using Type = float;
    using type = float;
};

template<>
struct BuildToValueMap<Fp16>
{
    using Type = float;
    using type = float;
};

template<>
struct BuildToValueMap<FpN>
{
    using Type = float;
    using type = float;
};

template<>
struct BuildToValueMap<Point>
{
    using Type = uint64_t;
    using type = uint64_t;
};

// --------------------------> utility functions related to alignment <------------------------------------

/// @brief return true if the specified pointer is aligned
__hostdev__ inline static bool isAligned(const void* p)
{
    return uint64_t(p) % NANOVDB_DATA_ALIGNMENT == 0;
}

/// @brief return true if the specified pointer is aligned and not NULL
__hostdev__ inline static bool isValid(const void* p)
{
    return p != nullptr && uint64_t(p) % NANOVDB_DATA_ALIGNMENT == 0;
}

/// @brief return the smallest number of bytes that when added to the specified pointer results in an aligned pointer
__hostdev__ inline static uint64_t alignmentPadding(const void* p)
{
    NANOVDB_ASSERT(p);
    return (NANOVDB_DATA_ALIGNMENT - (uint64_t(p) % NANOVDB_DATA_ALIGNMENT)) % NANOVDB_DATA_ALIGNMENT;
}

/// @brief offset the specified pointer so it is aligned.
template <typename T>
__hostdev__ inline static T* alignPtr(T* p)
{
    NANOVDB_ASSERT(p);
    return reinterpret_cast<T*>( (uint8_t*)p + alignmentPadding(p) );
}

/// @brief offset the specified const pointer so it is aligned.
template <typename T>
__hostdev__ inline static const T* alignPtr(const T* p)
{
    NANOVDB_ASSERT(p);
    return reinterpret_cast<const T*>( (const uint8_t*)p + alignmentPadding(p) );
}

// --------------------------> PtrDiff <------------------------------------

/// @brief Compute the distance, in bytes, between two pointers
/// @tparam T1 Type of the first pointer
/// @tparam T2 Type of the second pointer
/// @param p fist pointer, assumed to NOT be NULL
/// @param q second pointer, assumed to NOT be NULL
/// @return signed distance between pointer addresses in units of bytes
template<typename T1, typename T2>
__hostdev__ inline static int64_t PtrDiff(const T1* p, const T2* q)
{
    NANOVDB_ASSERT(p && q);
    return reinterpret_cast<const char*>(p) - reinterpret_cast<const char*>(q);
}

// --------------------------> PtrAdd <------------------------------------

/// @brief Adds a byte offset of a non-const pointer to produce another non-const pointer
/// @tparam DstT Type of the return pointer
/// @tparam SrcT Type of the input pointer
/// @param p non-const input pointer, assumed to NOT be NULL
/// @param offset signed byte offset
/// @return a non-const pointer defined as the offset of an input pointer
template<typename DstT, typename SrcT>
__hostdev__ inline static DstT* PtrAdd(SrcT* p, int64_t offset)
{
    NANOVDB_ASSERT(p);
    return reinterpret_cast<DstT*>(reinterpret_cast<char*>(p) + offset);
}

/// @brief Adds a byte offset of a const pointer to produce another const pointer
/// @tparam DstT Type of the return pointer
/// @tparam SrcT Type of the input pointer
/// @param p const input pointer, assumed to NOT be NULL
/// @param offset signed byte offset
/// @return a const pointer defined as the offset of a const input pointer
template<typename DstT, typename SrcT>
__hostdev__ inline static const DstT* PtrAdd(const SrcT* p, int64_t offset)
{
    NANOVDB_ASSERT(p);
    return reinterpret_cast<const DstT*>(reinterpret_cast<const char*>(p) + offset);
}

// --------------------------> isFloatingPoint(GridType) <------------------------------------

/// @brief return true if the GridType maps to a floating point type
__hostdev__ inline bool isFloatingPoint(GridType gridType)
{
    return gridType == GridType::Float ||
           gridType == GridType::Double ||
           gridType == GridType::Half ||
           gridType == GridType::Fp4 ||
           gridType == GridType::Fp8 ||
           gridType == GridType::Fp16 ||
           gridType == GridType::FpN;
}

// --------------------------> isFloatingPointVector(GridType) <------------------------------------

/// @brief return true if the GridType maps to a floating point vec3.
__hostdev__ inline bool isFloatingPointVector(GridType gridType)
{
    return gridType == GridType::Vec3f ||
           gridType == GridType::Vec3d ||
           gridType == GridType::Vec4f ||
           gridType == GridType::Vec4d;
}

// --------------------------> isInteger(GridType) <------------------------------------

/// @brief Return true if the GridType maps to a POD integer type.
/// @details These types are used to associate a voxel with a POD integer type
__hostdev__ inline bool isInteger(GridType gridType)
{
    return gridType == GridType::Int16 ||
           gridType == GridType::Int32 ||
           gridType == GridType::Int64 ||
           gridType == GridType::UInt32;
}

// --------------------------> isIndex(GridType) <------------------------------------

/// @brief Return true if the GridType maps to a special index type (not a POD integer type).
/// @details These types are used to index from a voxel into an external array of values, e.g. sidecar or blind data.
__hostdev__ inline bool isIndex(GridType gridType)
{
    return gridType == GridType::Index ||// index both active and inactive values
           gridType == GridType::OnIndex ||// index active values only
           gridType == GridType::IndexMask ||// as Index, but with an additional mask
           gridType == GridType::OnIndexMask;// as OnIndex, but with an additional mask
}

// --------------------------> memcpy64 <------------------------------------

/// @brief copy 64 bit words from @c src to @c dst
/// @param dst 64 bit aligned pointer to destination
/// @param src 64 bit aligned pointer to source
/// @param word_count number of 64 bit words to be copied
/// @return destination pointer @c dst
/// @warning @c src and @c dst cannot overlap and should both be 64 bit aligned
__hostdev__ inline static void* memcpy64(void *dst, const void *src, size_t word_count)
{
    NANOVDB_ASSERT(uint64_t(dst) % 8 == 0 && uint64_t(src) % 8 == 0);
    auto *d = reinterpret_cast<uint64_t*>(dst), *e = d + word_count;
    auto *s = reinterpret_cast<const uint64_t*>(src);
    while (d != e) *d++ = *s++;
    return dst;
}

// --------------------------> isValue(GridType, GridClass) <------------------------------------

/// @brief return true if the combination of GridType and GridClass is valid.
__hostdev__ inline bool isValid(GridType gridType, GridClass gridClass)
{
    if (gridClass == GridClass::LevelSet || gridClass == GridClass::FogVolume) {
        return isFloatingPoint(gridType);
    } else if (gridClass == GridClass::Staggered) {
        return isFloatingPointVector(gridType);
    } else if (gridClass == GridClass::PointIndex || gridClass == GridClass::PointData) {
        return gridType == GridType::PointIndex || gridType == GridType::UInt32;
    } else if (gridClass == GridClass::Topology) {
        return gridType == GridType::Mask;
    } else if (gridClass == GridClass::IndexGrid) {
        return isIndex(gridType);
    } else if (gridClass == GridClass::VoxelVolume) {
        return gridType == GridType::RGBA8 || gridType == GridType::Float ||
               gridType == GridType::Double || gridType == GridType::Vec3f ||
               gridType == GridType::Vec3d || gridType == GridType::UInt32;
    }
    return gridClass < GridClass::End && gridType < GridType::End; // any valid combination
}

// --------------------------> validation of blind data meta data <------------------------------------

/// @brief return true if the combination of GridBlindDataClass, GridBlindDataSemantic and GridType is valid.
__hostdev__ inline bool isValid(const GridBlindDataClass&    blindClass,
                                const GridBlindDataSemantic& blindSemantics,
                                const GridType&              blindType)
{
    bool test = false;
    switch (blindClass) {
    case GridBlindDataClass::IndexArray:
        test = (blindSemantics == GridBlindDataSemantic::Unknown ||
                blindSemantics == GridBlindDataSemantic::PointId) &&
               isInteger(blindType);
        break;
    case GridBlindDataClass::AttributeArray:
        if (blindSemantics == GridBlindDataSemantic::PointPosition ||
            blindSemantics == GridBlindDataSemantic::WorldCoords) {
            test = blindType == GridType::Vec3f || blindType == GridType::Vec3d;
        } else if (blindSemantics == GridBlindDataSemantic::GridCoords) {
            test = blindType == GridType::Vec3f;
        } else if (blindSemantics == GridBlindDataSemantic::VoxelCoords) {
            test = blindType == GridType::Vec3f || blindType == GridType::Vec3u8 || blindType == GridType::Vec3u16;
        } else {
            test = blindSemantics != GridBlindDataSemantic::PointId;
        }
        break;
    case GridBlindDataClass::GridName:
        test = blindSemantics == GridBlindDataSemantic::Unknown && blindType == GridType::Unknown;
        break;
    default: // captures blindClass == Unknown and ChannelArray
        test = blindClass < GridBlindDataClass::End &&
               blindSemantics < GridBlindDataSemantic::End &&
               blindType < GridType::End; // any valid combination
        break;
    }
    //if (!test) printf("Invalid combination: GridBlindDataClass=%u, GridBlindDataSemantic=%u, GridType=%u\n",(uint32_t)blindClass, (uint32_t)blindSemantics, (uint32_t)blindType);
    return test;
}

// ----------------------------> Version class <-------------------------------------

/// @brief Bit-compacted representation of all three version numbers
///
/// @details major is the top 11 bits, minor is the 11 middle bits and patch is the lower 10 bits
class Version
{
    uint32_t mData; // 11 + 11 + 10 bit packing of major + minor + patch
public:
    /// @brief Default constructor
    __hostdev__ Version()
        : mData(uint32_t(NANOVDB_MAJOR_VERSION_NUMBER) << 21 |
                uint32_t(NANOVDB_MINOR_VERSION_NUMBER) << 10 |
                uint32_t(NANOVDB_PATCH_VERSION_NUMBER))
    {
    }
    /// @brief Constructor from a raw uint32_t data representation
    __hostdev__ Version(uint32_t data) : mData(data) {}
    /// @brief Constructor from major.minor.patch version numbers
    __hostdev__ Version(uint32_t major, uint32_t minor, uint32_t patch)
        : mData(major << 21 | minor << 10 | patch)
    {
        NANOVDB_ASSERT(major < (1u << 11)); // max value of major is 2047
        NANOVDB_ASSERT(minor < (1u << 11)); // max value of minor is 2047
        NANOVDB_ASSERT(patch < (1u << 10)); // max value of patch is 1023
    }
    __hostdev__ bool     operator==(const Version& rhs) const { return mData == rhs.mData; }
    __hostdev__ bool     operator<( const Version& rhs) const { return mData < rhs.mData; }
    __hostdev__ bool     operator<=(const Version& rhs) const { return mData <= rhs.mData; }
    __hostdev__ bool     operator>( const Version& rhs) const { return mData > rhs.mData; }
    __hostdev__ bool     operator>=(const Version& rhs) const { return mData >= rhs.mData; }
    __hostdev__ uint32_t id() const { return mData; }
    __hostdev__ uint32_t getMajor() const { return (mData >> 21) & ((1u << 11) - 1); }
    __hostdev__ uint32_t getMinor() const { return (mData >> 10) & ((1u << 11) - 1); }
    __hostdev__ uint32_t getPatch() const { return  mData        & ((1u << 10) - 1); }
    __hostdev__ bool isCompatible() const { return this->getMajor() == uint32_t(NANOVDB_MAJOR_VERSION_NUMBER); }
    /// @brief Returns the difference between major version of this instance and NANOVDB_MAJOR_VERSION_NUMBER
    /// @return return 0 if the major version equals NANOVDB_MAJOR_VERSION_NUMBER, else a negative age if this
    ///         instance has a smaller major verion (is older), and a positive age if it is newer, i.e. larger.
    __hostdev__ int age() const {return int(this->getMajor()) - int(NANOVDB_MAJOR_VERSION_NUMBER);}

#ifndef __CUDACC_RTC__
    /// @brief returns a c-string of the semantic version, i.e. major.minor.patch
    const char* c_str() const
    {
        char* buffer = (char*)malloc(4 + 1 + 4 + 1 + 4 + 1); // xxxx.xxxx.xxxx\0
        snprintf(buffer, 4 + 1 + 4 + 1 + 4 + 1, "%u.%u.%u", this->getMajor(), this->getMinor(), this->getPatch()); // Prevents overflows by enforcing a fixed size of buffer
        return buffer;
    }
#endif
}; // Version

// ----------------------------> Various math functions <-------------------------------------

//@{
/// @brief Pi constant taken from Boost to match old behaviour
template<typename T>
inline __hostdev__ constexpr T pi()
{
    return 3.141592653589793238462643383279502884e+00;
}
template<>
inline __hostdev__ constexpr float pi()
{
    return 3.141592653589793238462643383279502884e+00F;
}
template<>
inline __hostdev__ constexpr double pi()
{
    return 3.141592653589793238462643383279502884e+00;
}
template<>
inline __hostdev__ constexpr long double pi()
{
    return 3.141592653589793238462643383279502884e+00L;
}
//@}

//@{
/// Tolerance for floating-point comparison
template<typename T>
struct Tolerance;
template<>
struct Tolerance<float>
{
    __hostdev__ static float value() { return 1e-8f; }
};
template<>
struct Tolerance<double>
{
    __hostdev__ static double value() { return 1e-15; }
};
//@}

//@{
/// Delta for small floating-point offsets
template<typename T>
struct Delta;
template<>
struct Delta<float>
{
    __hostdev__ static float value() { return 1e-5f; }
};
template<>
struct Delta<double>
{
    __hostdev__ static double value() { return 1e-9; }
};
//@}

//@{
/// Maximum floating-point values
template<typename T>
struct Maximum;
#if defined(__CUDA_ARCH__) || defined(__HIP__)
template<>
struct Maximum<int>
{
    __hostdev__ static int value() { return 2147483647; }
};
template<>
struct Maximum<uint32_t>
{
    __hostdev__ static uint32_t value() { return 4294967295u; }
};
template<>
struct Maximum<float>
{
    __hostdev__ static float value() { return 1e+38f; }
};
template<>
struct Maximum<double>
{
    __hostdev__ static double value() { return 1e+308; }
};
#else
template<typename T>
struct Maximum
{
    static T value() { return std::numeric_limits<T>::max(); }
};
#endif
//@}

template<typename Type>
__hostdev__ inline bool isApproxZero(const Type& x)
{
    return !(x > Tolerance<Type>::value()) && !(x < -Tolerance<Type>::value());
}

template<typename Type>
__hostdev__ inline Type Min(Type a, Type b)
{
    return (a < b) ? a : b;
}
__hostdev__ inline int32_t Min(int32_t a, int32_t b)
{
    return int32_t(fminf(float(a), float(b)));
}
__hostdev__ inline uint32_t Min(uint32_t a, uint32_t b)
{
    return uint32_t(fminf(float(a), float(b)));
}
__hostdev__ inline float Min(float a, float b)
{
    return fminf(a, b);
}
__hostdev__ inline double Min(double a, double b)
{
    return fmin(a, b);
}
template<typename Type>
__hostdev__ inline Type Max(Type a, Type b)
{
    return (a > b) ? a : b;
}

__hostdev__ inline int32_t Max(int32_t a, int32_t b)
{
    return int32_t(fmaxf(float(a), float(b)));
}
__hostdev__ inline uint32_t Max(uint32_t a, uint32_t b)
{
    return uint32_t(fmaxf(float(a), float(b)));
}
__hostdev__ inline float Max(float a, float b)
{
    return fmaxf(a, b);
}
__hostdev__ inline double Max(double a, double b)
{
    return fmax(a, b);
}
__hostdev__ inline float Clamp(float x, float a, float b)
{
    return Max(Min(x, b), a);
}
__hostdev__ inline double Clamp(double x, double a, double b)
{
    return Max(Min(x, b), a);
}

__hostdev__ inline float Fract(float x)
{
    return x - floorf(x);
}
__hostdev__ inline double Fract(double x)
{
    return x - floor(x);
}

__hostdev__ inline int32_t Floor(float x)
{
    return int32_t(floorf(x));
}
__hostdev__ inline int32_t Floor(double x)
{
    return int32_t(floor(x));
}

__hostdev__ inline int32_t Ceil(float x)
{
    return int32_t(ceilf(x));
}
__hostdev__ inline int32_t Ceil(double x)
{
    return int32_t(ceil(x));
}

template<typename T>
__hostdev__ inline T Pow2(T x)
{
    return x * x;
}

template<typename T>
__hostdev__ inline T Pow3(T x)
{
    return x * x * x;
}

template<typename T>
__hostdev__ inline T Pow4(T x)
{
    return Pow2(x * x);
}
template<typename T>
__hostdev__ inline T Abs(T x)
{
    return x < 0 ? -x : x;
}

template<>
__hostdev__ inline float Abs(float x)
{
    return fabsf(x);
}

template<>
__hostdev__ inline double Abs(double x)
{
    return fabs(x);
}

template<>
__hostdev__ inline int Abs(int x)
{
    return abs(x);
}

template<typename CoordT, typename RealT, template<typename> class Vec3T>
__hostdev__ inline CoordT Round(const Vec3T<RealT>& xyz);

template<typename CoordT, template<typename> class Vec3T>
__hostdev__ inline CoordT Round(const Vec3T<float>& xyz)
{
    return CoordT(int32_t(rintf(xyz[0])), int32_t(rintf(xyz[1])), int32_t(rintf(xyz[2])));
    //return CoordT(int32_t(roundf(xyz[0])), int32_t(roundf(xyz[1])), int32_t(roundf(xyz[2])) );
    //return CoordT(int32_t(floorf(xyz[0] + 0.5f)), int32_t(floorf(xyz[1] + 0.5f)), int32_t(floorf(xyz[2] + 0.5f)));
}

template<typename CoordT, template<typename> class Vec3T>
__hostdev__ inline CoordT Round(const Vec3T<double>& xyz)
{
    return CoordT(int32_t(floor(xyz[0] + 0.5)), int32_t(floor(xyz[1] + 0.5)), int32_t(floor(xyz[2] + 0.5)));
}

template<typename CoordT, typename RealT, template<typename> class Vec3T>
__hostdev__ inline CoordT RoundDown(const Vec3T<RealT>& xyz)
{
    return CoordT(Floor(xyz[0]), Floor(xyz[1]), Floor(xyz[2]));
}

//@{
/// Return the square root of a floating-point value.
__hostdev__ inline float Sqrt(float x)
{
    return sqrtf(x);
}
__hostdev__ inline double Sqrt(double x)
{
    return sqrt(x);
}
//@}

/// Return the sign of the given value as an integer (either -1, 0 or 1).
template<typename T>
__hostdev__ inline T Sign(const T& x)
{
    return ((T(0) < x) ? T(1) : T(0)) - ((x < T(0)) ? T(1) : T(0));
}

template<typename Vec3T>
__hostdev__ inline int MinIndex(const Vec3T& v)
{
#if 0
    static const int hashTable[8] = {2, 1, 9, 1, 2, 9, 0, 0}; //9 are dummy values
    const int        hashKey = ((v[0] < v[1]) << 2) + ((v[0] < v[2]) << 1) + (v[1] < v[2]); // ?*4+?*2+?*1
    return hashTable[hashKey];
#else
    if (v[0] < v[1] && v[0] < v[2])
        return 0;
    if (v[1] < v[2])
        return 1;
    else
        return 2;
#endif
}

template<typename Vec3T>
__hostdev__ inline int MaxIndex(const Vec3T& v)
{
#if 0
    static const int hashTable[8] = {2, 1, 9, 1, 2, 9, 0, 0}; //9 are dummy values
    const int        hashKey = ((v[0] > v[1]) << 2) + ((v[0] > v[2]) << 1) + (v[1] > v[2]); // ?*4+?*2+?*1
    return hashTable[hashKey];
#else
    if (v[0] > v[1] && v[0] > v[2])
        return 0;
    if (v[1] > v[2])
        return 1;
    else
        return 2;
#endif
}

/// @brief round up byteSize to the nearest wordSize, e.g. to align to machine word: AlignUp<sizeof(size_t)(n)
///
/// @details both wordSize and byteSize are in byte units
template<uint64_t wordSize>
__hostdev__ inline uint64_t AlignUp(uint64_t byteCount)
{
    const uint64_t r = byteCount % wordSize;
    return r ? byteCount - r + wordSize : byteCount;
}

// ------------------------------> Coord <--------------------------------------

// forward declaration so we can define Coord::asVec3s and Coord::asVec3d
template<typename>
class Vec3;

/// @brief Signed (i, j, k) 32-bit integer coordinate class, similar to openvdb::math::Coord
class Coord
{
    int32_t mVec[3]; // private member data - three signed index coordinates
public:
    using ValueType = int32_t;
    using IndexType = uint32_t;

    /// @brief Initialize all coordinates to zero.
    __hostdev__ Coord()
        : mVec{0, 0, 0}
    {
    }

    /// @brief Initializes all coordinates to the given signed integer.
    __hostdev__ explicit Coord(ValueType n)
        : mVec{n, n, n}
    {
    }

    /// @brief Initializes coordinate to the given signed integers.
    __hostdev__ Coord(ValueType i, ValueType j, ValueType k)
        : mVec{i, j, k}
    {
    }

    __hostdev__ Coord(ValueType* ptr)
        : mVec{ptr[0], ptr[1], ptr[2]}
    {
    }

    __hostdev__ int32_t x() const { return mVec[0]; }
    __hostdev__ int32_t y() const { return mVec[1]; }
    __hostdev__ int32_t z() const { return mVec[2]; }

    __hostdev__ int32_t& x() { return mVec[0]; }
    __hostdev__ int32_t& y() { return mVec[1]; }
    __hostdev__ int32_t& z() { return mVec[2]; }

    __hostdev__ static Coord max() { return Coord(int32_t((1u << 31) - 1)); }

    __hostdev__ static Coord min() { return Coord(-int32_t((1u << 31) - 1) - 1); }

    __hostdev__ static size_t memUsage() { return sizeof(Coord); }

    /// @brief Return a const reference to the given Coord component.
    /// @warning The argument is assumed to be 0, 1, or 2.
    __hostdev__ const ValueType& operator[](IndexType i) const { return mVec[i]; }

    /// @brief Return a non-const reference to the given Coord component.
    /// @warning The argument is assumed to be 0, 1, or 2.
    __hostdev__ ValueType& operator[](IndexType i) { return mVec[i]; }

    /// @brief Assignment operator that works with openvdb::Coord
    template<typename CoordT>
    __hostdev__ Coord& operator=(const CoordT& other)
    {
        static_assert(sizeof(Coord) == sizeof(CoordT), "Mis-matched sizeof");
        mVec[0] = other[0];
        mVec[1] = other[1];
        mVec[2] = other[2];
        return *this;
    }

    /// @brief Return a new instance with coordinates masked by the given unsigned integer.
    __hostdev__ Coord operator&(IndexType n) const { return Coord(mVec[0] & n, mVec[1] & n, mVec[2] & n); }

    // @brief Return a new instance with coordinates left-shifted by the given unsigned integer.
    __hostdev__ Coord operator<<(IndexType n) const { return Coord(mVec[0] << n, mVec[1] << n, mVec[2] << n); }

    // @brief Return a new instance with coordinates right-shifted by the given unsigned integer.
    __hostdev__ Coord operator>>(IndexType n) const { return Coord(mVec[0] >> n, mVec[1] >> n, mVec[2] >> n); }

    /// @brief Return true if this Coord is lexicographically less than the given Coord.
    __hostdev__ bool operator<(const Coord& rhs) const
    {
        return mVec[0] < rhs[0] ? true
             : mVec[0] > rhs[0] ? false
             : mVec[1] < rhs[1] ? true
             : mVec[1] > rhs[1] ? false
             : mVec[2] < rhs[2] ? true : false;
    }

    /// @brief Return true if this Coord is lexicographically less or equal to the given Coord.
    __hostdev__ bool operator<=(const Coord& rhs) const
    {
        return mVec[0] < rhs[0] ? true
             : mVec[0] > rhs[0] ? false
             : mVec[1] < rhs[1] ? true
             : mVec[1] > rhs[1] ? false
             : mVec[2] <=rhs[2] ? true : false;
    }

    // @brief Return true if the Coord components are identical.
    __hostdev__ bool   operator==(const Coord& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; }
    __hostdev__ bool   operator!=(const Coord& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; }
    __hostdev__ Coord& operator&=(int n)
    {
        mVec[0] &= n;
        mVec[1] &= n;
        mVec[2] &= n;
        return *this;
    }
    __hostdev__ Coord& operator<<=(uint32_t n)
    {
        mVec[0] <<= n;
        mVec[1] <<= n;
        mVec[2] <<= n;
        return *this;
    }
    __hostdev__ Coord& operator>>=(uint32_t n)
    {
        mVec[0] >>= n;
        mVec[1] >>= n;
        mVec[2] >>= n;
        return *this;
    }
    __hostdev__ Coord& operator+=(int n)
    {
        mVec[0] += n;
        mVec[1] += n;
        mVec[2] += n;
        return *this;
    }
    __hostdev__ Coord  operator+(const Coord& rhs) const { return Coord(mVec[0] + rhs[0], mVec[1] + rhs[1], mVec[2] + rhs[2]); }
    __hostdev__ Coord  operator-(const Coord& rhs) const { return Coord(mVec[0] - rhs[0], mVec[1] - rhs[1], mVec[2] - rhs[2]); }
    __hostdev__ Coord  operator-() const { return Coord(-mVec[0], -mVec[1], -mVec[2]); }
    __hostdev__ Coord& operator+=(const Coord& rhs)
    {
        mVec[0] += rhs[0];
        mVec[1] += rhs[1];
        mVec[2] += rhs[2];
        return *this;
    }
    __hostdev__ Coord& operator-=(const Coord& rhs)
    {
        mVec[0] -= rhs[0];
        mVec[1] -= rhs[1];
        mVec[2] -= rhs[2];
        return *this;
    }

    /// @brief Perform a component-wise minimum with the other Coord.
    __hostdev__ Coord& minComponent(const Coord& other)
    {
        if (other[0] < mVec[0])
            mVec[0] = other[0];
        if (other[1] < mVec[1])
            mVec[1] = other[1];
        if (other[2] < mVec[2])
            mVec[2] = other[2];
        return *this;
    }

    /// @brief Perform a component-wise maximum with the other Coord.
    __hostdev__ Coord& maxComponent(const Coord& other)
    {
        if (other[0] > mVec[0])
            mVec[0] = other[0];
        if (other[1] > mVec[1])
            mVec[1] = other[1];
        if (other[2] > mVec[2])
            mVec[2] = other[2];
        return *this;
    }
#if defined(__CUDACC__) // the following functions only run on the GPU!
    __device__ inline Coord& minComponentAtomic(const Coord& other)
    {
        atomicMin(&mVec[0], other[0]);
        atomicMin(&mVec[1], other[1]);
        atomicMin(&mVec[2], other[2]);
        return *this;
    }
    __device__ inline Coord& maxComponentAtomic(const Coord& other)
    {
        atomicMax(&mVec[0], other[0]);
        atomicMax(&mVec[1], other[1]);
        atomicMax(&mVec[2], other[2]);
        return *this;
    }
#endif

    __hostdev__ Coord offsetBy(ValueType dx, ValueType dy, ValueType dz) const
    {
        return Coord(mVec[0] + dx, mVec[1] + dy, mVec[2] + dz);
    }

    __hostdev__ Coord offsetBy(ValueType n) const { return this->offsetBy(n, n, n); }

    /// Return true if any of the components of @a a are smaller than the
    /// corresponding components of @a b.
    __hostdev__ static inline bool lessThan(const Coord& a, const Coord& b)
    {
        return (a[0] < b[0] || a[1] < b[1] || a[2] < b[2]);
    }

    /// @brief Return the largest integer coordinates that are not greater
    /// than @a xyz (node centered conversion).
    template<typename Vec3T>
    __hostdev__ static Coord Floor(const Vec3T& xyz) { return Coord(nanovdb::Floor(xyz[0]), nanovdb::Floor(xyz[1]), nanovdb::Floor(xyz[2])); }

    /// @brief Return a hash key derived from the existing coordinates.
    /// @details The hash function is originally taken from the SIGGRAPH paper:
    ///          "VDB: High-resolution sparse volumes with dynamic topology"
    ///          and the prime numbers are modified based on the ACM Transactions on Graphics paper:
    ///          "Real-time 3D reconstruction at scale using voxel hashing" (the second number had a typo!)
    template<int Log2N = 3 + 4 + 5>
    __hostdev__ uint32_t hash() const { return ((1 << Log2N) - 1) & (mVec[0] * 73856093 ^ mVec[1] * 19349669 ^ mVec[2] * 83492791); }

    /// @brief Return the octant of this Coord
    //__hostdev__ size_t octant() const { return (uint32_t(mVec[0])>>31) | ((uint32_t(mVec[1])>>31)<<1) | ((uint32_t(mVec[2])>>31)<<2); }
    __hostdev__ uint8_t octant() const { return (uint8_t(bool(mVec[0] & (1u << 31)))) |
                                                (uint8_t(bool(mVec[1] & (1u << 31))) << 1) |
                                                (uint8_t(bool(mVec[2] & (1u << 31))) << 2); }

    /// @brief Return a single precision floating-point vector of this coordinate
    __hostdev__ inline Vec3<float> asVec3s() const;

    /// @brief Return a double precision floating-point vector of this coordinate
    __hostdev__ inline Vec3<double> asVec3d() const;

    // returns a copy of itself, so it mimics the behaviour of Vec3<T>::round()
    __hostdev__ inline Coord round() const { return *this; }
}; // Coord class

// ----------------------------> Vec3 <--------------------------------------

/// @brief A simple vector class with three components, similar to openvdb::math::Vec3
template<typename T>
class Vec3
{
    T mVec[3];

public:
    static const int SIZE = 3;
    static const int size = 3; // in openvdb::math::Tuple
    using ValueType = T;
    Vec3() = default;
    __hostdev__ explicit Vec3(T x)
        : mVec{x, x, x}
    {
    }
    __hostdev__ Vec3(T x, T y, T z)
        : mVec{x, y, z}
    {
    }
    template<template<class> class Vec3T, class T2>
    __hostdev__ Vec3(const Vec3T<T2>& v)
        : mVec{T(v[0]), T(v[1]), T(v[2])}
    {
        static_assert(Vec3T<T2>::size == size, "expected Vec3T::size==3!");
    }
    template<typename T2>
    __hostdev__ explicit Vec3(const Vec3<T2>& v)
        : mVec{T(v[0]), T(v[1]), T(v[2])}
    {
    }
    __hostdev__ explicit Vec3(const Coord& ijk)
        : mVec{T(ijk[0]), T(ijk[1]), T(ijk[2])}
    {
    }
    __hostdev__ bool operator==(const Vec3& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; }
    __hostdev__ bool operator!=(const Vec3& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; }
    template<template<class> class Vec3T, class T2>
    __hostdev__ Vec3& operator=(const Vec3T<T2>& rhs)
    {
        static_assert(Vec3T<T2>::size == size, "expected Vec3T::size==3!");
        mVec[0] = rhs[0];
        mVec[1] = rhs[1];
        mVec[2] = rhs[2];
        return *this;
    }
    __hostdev__ const T& operator[](int i) const { return mVec[i]; }
    __hostdev__ T&       operator[](int i) { return mVec[i]; }
    template<typename Vec3T>
    __hostdev__ T dot(const Vec3T& v) const { return mVec[0] * v[0] + mVec[1] * v[1] + mVec[2] * v[2]; }
    template<typename Vec3T>
    __hostdev__ Vec3 cross(const Vec3T& v) const
    {
        return Vec3(mVec[1] * v[2] - mVec[2] * v[1],
                    mVec[2] * v[0] - mVec[0] * v[2],
                    mVec[0] * v[1] - mVec[1] * v[0]);
    }
    __hostdev__ T lengthSqr() const
    {
        return mVec[0] * mVec[0] + mVec[1] * mVec[1] + mVec[2] * mVec[2]; // 5 flops
    }
    __hostdev__ T     length() const { return Sqrt(this->lengthSqr()); }
    __hostdev__ Vec3  operator-() const { return Vec3(-mVec[0], -mVec[1], -mVec[2]); }
    __hostdev__ Vec3  operator*(const Vec3& v) const { return Vec3(mVec[0] * v[0], mVec[1] * v[1], mVec[2] * v[2]); }
    __hostdev__ Vec3  operator/(const Vec3& v) const { return Vec3(mVec[0] / v[0], mVec[1] / v[1], mVec[2] / v[2]); }
    __hostdev__ Vec3  operator+(const Vec3& v) const { return Vec3(mVec[0] + v[0], mVec[1] + v[1], mVec[2] + v[2]); }
    __hostdev__ Vec3  operator-(const Vec3& v) const { return Vec3(mVec[0] - v[0], mVec[1] - v[1], mVec[2] - v[2]); }
    __hostdev__ Vec3  operator+(const Coord& ijk) const { return Vec3(mVec[0] + ijk[0], mVec[1] + ijk[1], mVec[2] + ijk[2]); }
    __hostdev__ Vec3  operator-(const Coord& ijk) const { return Vec3(mVec[0] - ijk[0], mVec[1] - ijk[1], mVec[2] - ijk[2]); }
    __hostdev__ Vec3  operator*(const T& s) const { return Vec3(s * mVec[0], s * mVec[1], s * mVec[2]); }
    __hostdev__ Vec3  operator/(const T& s) const { return (T(1) / s) * (*this); }
    __hostdev__ Vec3& operator+=(const Vec3& v)
    {
        mVec[0] += v[0];
        mVec[1] += v[1];
        mVec[2] += v[2];
        return *this;
    }
    __hostdev__ Vec3& operator+=(const Coord& ijk)
    {
        mVec[0] += T(ijk[0]);
        mVec[1] += T(ijk[1]);
        mVec[2] += T(ijk[2]);
        return *this;
    }
    __hostdev__ Vec3& operator-=(const Vec3& v)
    {
        mVec[0] -= v[0];
        mVec[1] -= v[1];
        mVec[2] -= v[2];
        return *this;
    }
    __hostdev__ Vec3& operator-=(const Coord& ijk)
    {
        mVec[0] -= T(ijk[0]);
        mVec[1] -= T(ijk[1]);
        mVec[2] -= T(ijk[2]);
        return *this;
    }
    __hostdev__ Vec3& operator*=(const T& s)
    {
        mVec[0] *= s;
        mVec[1] *= s;
        mVec[2] *= s;
        return *this;
    }
    __hostdev__ Vec3& operator/=(const T& s) { return (*this) *= T(1) / s; }
    __hostdev__ Vec3& normalize() { return (*this) /= this->length(); }
    /// @brief Perform a component-wise minimum with the other Coord.
    __hostdev__ Vec3& minComponent(const Vec3& other)
    {
        if (other[0] < mVec[0])
            mVec[0] = other[0];
        if (other[1] < mVec[1])
            mVec[1] = other[1];
        if (other[2] < mVec[2])
            mVec[2] = other[2];
        return *this;
    }

    /// @brief Perform a component-wise maximum with the other Coord.
    __hostdev__ Vec3& maxComponent(const Vec3& other)
    {
        if (other[0] > mVec[0])
            mVec[0] = other[0];
        if (other[1] > mVec[1])
            mVec[1] = other[1];
        if (other[2] > mVec[2])
            mVec[2] = other[2];
        return *this;
    }
    /// @brief Return the smallest vector component
    __hostdev__ ValueType min() const
    {
        return mVec[0] < mVec[1] ? (mVec[0] < mVec[2] ? mVec[0] : mVec[2]) : (mVec[1] < mVec[2] ? mVec[1] : mVec[2]);
    }
    /// @brief Return the largest vector component
    __hostdev__ ValueType max() const
    {
        return mVec[0] > mVec[1] ? (mVec[0] > mVec[2] ? mVec[0] : mVec[2]) : (mVec[1] > mVec[2] ? mVec[1] : mVec[2]);
    }
    /// @brief Round each component if this Vec<T> up to its integer value
    /// @return Return an integer Coord
    __hostdev__ Coord floor() const { return Coord(Floor(mVec[0]), Floor(mVec[1]), Floor(mVec[2])); }
    /// @brief Round each component if this Vec<T> down to its integer value
    /// @return Return an integer Coord
    __hostdev__ Coord ceil() const { return Coord(Ceil(mVec[0]), Ceil(mVec[1]), Ceil(mVec[2])); }
    /// @brief Round each component if this Vec<T> to its closest integer value
    /// @return Return an integer Coord
    __hostdev__ Coord round() const
    {
        if constexpr(is_same<T, float>::value) {
            return Coord(Floor(mVec[0] + 0.5f), Floor(mVec[1] + 0.5f), Floor(mVec[2] + 0.5f));
        } else if constexpr(is_same<T, int>::value) {
            return Coord(mVec[0], mVec[1], mVec[2]);
        } else {
            return Coord(Floor(mVec[0] + 0.5), Floor(mVec[1] + 0.5), Floor(mVec[2] + 0.5));
        }
    }

    /// @brief return a non-const raw constant pointer to array of three vector components
    __hostdev__ T* asPointer() { return mVec; }
    /// @brief return a const raw constant pointer to array of three vector components
    __hostdev__ const T* asPointer() const { return mVec; }
}; // Vec3<T>

template<typename T1, typename T2>
__hostdev__ inline Vec3<T2> operator*(T1 scalar, const Vec3<T2>& vec)
{
    return Vec3<T2>(scalar * vec[0], scalar * vec[1], scalar * vec[2]);
}
template<typename T1, typename T2>
__hostdev__ inline Vec3<T2> operator/(T1 scalar, const Vec3<T2>& vec)
{
    return Vec3<T2>(scalar / vec[0], scalar / vec[1], scalar / vec[2]);
}

//using Vec3R = Vec3<double>;// deprecated
using Vec3d = Vec3<double>;
using Vec3f = Vec3<float>;
using Vec3i = Vec3<int32_t>;
using Vec3u = Vec3<uint32_t>;
using Vec3u8 = Vec3<uint8_t>;
using Vec3u16 = Vec3<uint16_t>;

/// @brief Return a single precision floating-point vector of this coordinate
__hostdev__ inline Vec3f Coord::asVec3s() const
{
    return Vec3f(float(mVec[0]), float(mVec[1]), float(mVec[2]));
}

/// @brief Return a double precision floating-point vector of this coordinate
__hostdev__ inline Vec3d Coord::asVec3d() const
{
    return Vec3d(double(mVec[0]), double(mVec[1]), double(mVec[2]));
}

// ----------------------------> Vec4 <--------------------------------------

/// @brief A simple vector class with four components, similar to openvdb::math::Vec4
template<typename T>
class Vec4
{
    T mVec[4];

public:
    static const int SIZE = 4;
    static const int size = 4;
    using ValueType = T;
    Vec4() = default;
    __hostdev__ explicit Vec4(T x)
        : mVec{x, x, x, x}
    {
    }
    __hostdev__ Vec4(T x, T y, T z, T w)
        : mVec{x, y, z, w}
    {
    }
    template<typename T2>
    __hostdev__ explicit Vec4(const Vec4<T2>& v)
        : mVec{T(v[0]), T(v[1]), T(v[2]), T(v[3])}
    {
    }
    template<template<class> class Vec4T, class T2>
    __hostdev__ Vec4(const Vec4T<T2>& v)
        : mVec{T(v[0]), T(v[1]), T(v[2]), T(v[3])}
    {
        static_assert(Vec4T<T2>::size == size, "expected Vec4T::size==4!");
    }
    __hostdev__ bool operator==(const Vec4& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2] && mVec[3] == rhs[3]; }
    __hostdev__ bool operator!=(const Vec4& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2] || mVec[3] != rhs[3]; }
    template<template<class> class Vec4T, class T2>
    __hostdev__ Vec4& operator=(const Vec4T<T2>& rhs)
    {
        static_assert(Vec4T<T2>::size == size, "expected Vec4T::size==4!");
        mVec[0] = rhs[0];
        mVec[1] = rhs[1];
        mVec[2] = rhs[2];
        mVec[3] = rhs[3];
        return *this;
    }

    __hostdev__ const T& operator[](int i) const { return mVec[i]; }
    __hostdev__ T&       operator[](int i) { return mVec[i]; }
    template<typename Vec4T>
    __hostdev__ T dot(const Vec4T& v) const { return mVec[0] * v[0] + mVec[1] * v[1] + mVec[2] * v[2] + mVec[3] * v[3]; }
    __hostdev__ T lengthSqr() const
    {
        return mVec[0] * mVec[0] + mVec[1] * mVec[1] + mVec[2] * mVec[2] + mVec[3] * mVec[3]; // 7 flops
    }
    __hostdev__ T     length() const { return Sqrt(this->lengthSqr()); }
    __hostdev__ Vec4  operator-() const { return Vec4(-mVec[0], -mVec[1], -mVec[2], -mVec[3]); }
    __hostdev__ Vec4  operator*(const Vec4& v) const { return Vec4(mVec[0] * v[0], mVec[1] * v[1], mVec[2] * v[2], mVec[3] * v[3]); }
    __hostdev__ Vec4  operator/(const Vec4& v) const { return Vec4(mVec[0] / v[0], mVec[1] / v[1], mVec[2] / v[2], mVec[3] / v[3]); }
    __hostdev__ Vec4  operator+(const Vec4& v) const { return Vec4(mVec[0] + v[0], mVec[1] + v[1], mVec[2] + v[2], mVec[3] + v[3]); }
    __hostdev__ Vec4  operator-(const Vec4& v) const { return Vec4(mVec[0] - v[0], mVec[1] - v[1], mVec[2] - v[2], mVec[3] - v[3]); }
    __hostdev__ Vec4  operator*(const T& s) const { return Vec4(s * mVec[0], s * mVec[1], s * mVec[2], s * mVec[3]); }
    __hostdev__ Vec4  operator/(const T& s) const { return (T(1) / s) * (*this); }
    __hostdev__ Vec4& operator+=(const Vec4& v)
    {
        mVec[0] += v[0];
        mVec[1] += v[1];
        mVec[2] += v[2];
        mVec[3] += v[3];
        return *this;
    }
    __hostdev__ Vec4& operator-=(const Vec4& v)
    {
        mVec[0] -= v[0];
        mVec[1] -= v[1];
        mVec[2] -= v[2];
        mVec[3] -= v[3];
        return *this;
    }
    __hostdev__ Vec4& operator*=(const T& s)
    {
        mVec[0] *= s;
        mVec[1] *= s;
        mVec[2] *= s;
        mVec[3] *= s;
        return *this;
    }
    __hostdev__ Vec4& operator/=(const T& s) { return (*this) *= T(1) / s; }
    __hostdev__ Vec4& normalize() { return (*this) /= this->length(); }
    /// @brief Perform a component-wise minimum with the other Coord.
    __hostdev__ Vec4& minComponent(const Vec4& other)
    {
        if (other[0] < mVec[0])
            mVec[0] = other[0];
        if (other[1] < mVec[1])
            mVec[1] = other[1];
        if (other[2] < mVec[2])
            mVec[2] = other[2];
        if (other[3] < mVec[3])
            mVec[3] = other[3];
        return *this;
    }

    /// @brief Perform a component-wise maximum with the other Coord.
    __hostdev__ Vec4& maxComponent(const Vec4& other)
    {
        if (other[0] > mVec[0])
            mVec[0] = other[0];
        if (other[1] > mVec[1])
            mVec[1] = other[1];
        if (other[2] > mVec[2])
            mVec[2] = other[2];
        if (other[3] > mVec[3])
            mVec[3] = other[3];
        return *this;
    }
}; // Vec4<T>

template<typename T1, typename T2>
__hostdev__ inline Vec4<T2> operator*(T1 scalar, const Vec4<T2>& vec)
{
    return Vec4<T2>(scalar * vec[0], scalar * vec[1], scalar * vec[2], scalar * vec[3]);
}
template<typename T1, typename T2>
__hostdev__ inline Vec4<T2> operator/(T1 scalar, const Vec4<T2>& vec)
{
    return Vec4<T2>(scalar / vec[0], scalar / vec[1], scalar / vec[2], scalar / vec[3]);
}

using Vec4R = Vec4<double>;
using Vec4d = Vec4<double>;
using Vec4f = Vec4<float>;
using Vec4i = Vec4<int>;


// --------------------------> Rgba8 <------------------------------------

/// @brief 8-bit red, green, blue, alpha packed into 32 bit unsigned int
class Rgba8
{
    union
    {
        uint8_t  c[4];   // 4 integer color channels of red, green, blue and alpha components.
        uint32_t packed; // 32 bit packed representation
    } mData;

public:
    static const int SIZE = 4;
    using ValueType = uint8_t;

    /// @brief Default copy constructor
    Rgba8(const Rgba8&) = default;

    /// @brief Default move constructor
    Rgba8(Rgba8&&) = default;

    /// @brief Default move assignment operator
    /// @return non-const reference to this instance
    Rgba8&      operator=(Rgba8&&) = default;

    /// @brief Default copy assignment operator
    /// @return non-const reference to this instance
    Rgba8&      operator=(const Rgba8&) = default;

    /// @brief Default ctor initializes all channels to zero
    __hostdev__ Rgba8()
        : mData{{0, 0, 0, 0}}
    {
        static_assert(sizeof(uint32_t) == sizeof(Rgba8), "Unexpected sizeof");
    }

    /// @brief integer r,g,b,a ctor where alpha channel defaults to opaque
    /// @note all values should be in the range 0u to 255u
    __hostdev__ Rgba8(uint8_t r, uint8_t g, uint8_t b, uint8_t a = 255u)
        : mData{{r, g, b, a}}
    {
    }

    /// @brief  @brief ctor where all channels are initialized to the same value
    /// @note value should be in the range 0u to 255u
    explicit __hostdev__ Rgba8(uint8_t v)
        : mData{{v, v, v, v}}
    {
    }

    /// @brief floating-point r,g,b,a ctor where alpha channel defaults to opaque
    /// @note all values should be in the range 0.0f to 1.0f
    __hostdev__ Rgba8(float r, float g, float b, float a = 1.0f)
        : mData{{static_cast<uint8_t>(0.5f + r * 255.0f), // round floats to nearest integers
                 static_cast<uint8_t>(0.5f + g * 255.0f), // double {{}} is needed due to union
                 static_cast<uint8_t>(0.5f + b * 255.0f),
                 static_cast<uint8_t>(0.5f + a * 255.0f)}}
    {
    }

    /// @brief Vec3f r,g,b ctor (alpha channel it set to 1)
    /// @note all values should be in the range 0.0f to 1.0f
    __hostdev__ Rgba8(const Vec3f& rgb)
        : Rgba8(rgb[0], rgb[1], rgb[2])
    {
    }

    /// @brief Vec4f r,g,b,a ctor
    /// @note all values should be in the range 0.0f to 1.0f
    __hostdev__ Rgba8(const Vec4f& rgba)
        : Rgba8(rgba[0], rgba[1], rgba[2], rgba[3])
    {
    }

    __hostdev__ bool  operator< (const Rgba8& rhs) const { return mData.packed < rhs.mData.packed; }
    __hostdev__ bool  operator==(const Rgba8& rhs) const { return mData.packed == rhs.mData.packed; }
    __hostdev__ float lengthSqr() const
    {
        return 0.0000153787005f * (float(mData.c[0]) * mData.c[0] +
                                   float(mData.c[1]) * mData.c[1] +
                                   float(mData.c[2]) * mData.c[2]); //1/255^2
    }
    __hostdev__ float           length() const { return sqrtf(this->lengthSqr()); }
    /// @brief return n'th color channel as a float in the range 0 to 1
    __hostdev__ float           asFloat(int n) const { return 0.003921569f*float(mData.c[n]); }// divide by 255
    __hostdev__ const uint8_t&  operator[](int n) const { return mData.c[n]; }
    __hostdev__ uint8_t&        operator[](int n) { return mData.c[n]; }
    __hostdev__ const uint32_t& packed() const { return mData.packed; }
    __hostdev__ uint32_t&       packed() { return mData.packed; }
    __hostdev__ const uint8_t&  r() const { return mData.c[0]; }
    __hostdev__ const uint8_t&  g() const { return mData.c[1]; }
    __hostdev__ const uint8_t&  b() const { return mData.c[2]; }
    __hostdev__ const uint8_t&  a() const { return mData.c[3]; }
    __hostdev__ uint8_t&        r() { return mData.c[0]; }
    __hostdev__ uint8_t&        g() { return mData.c[1]; }
    __hostdev__ uint8_t&        b() { return mData.c[2]; }
    __hostdev__ uint8_t&        a() { return mData.c[3]; }
    __hostdev__                 operator Vec3f() const {
        return Vec3f(this->asFloat(0), this->asFloat(1), this->asFloat(2));
    }
    __hostdev__                 operator Vec4f() const {
        return Vec4f(this->asFloat(0), this->asFloat(1), this->asFloat(2), this->asFloat(3));
    }
}; // Rgba8

using PackedRGBA8 = Rgba8; // for backwards compatibility

// ----------------------------> TensorTraits <--------------------------------------

template<typename T, int Rank = (is_specialization<T, Vec3>::value || is_specialization<T, Vec4>::value || is_same<T, Rgba8>::value) ? 1 : 0>
struct TensorTraits;

template<typename T>
struct TensorTraits<T, 0>
{
    static const int  Rank = 0; // i.e. scalar
    static const bool IsScalar = true;
    static const bool IsVector = false;
    static const int  Size = 1;
    using ElementType = T;
    static T scalar(const T& s) { return s; }
};

template<typename T>
struct TensorTraits<T, 1>
{
    static const int  Rank = 1; // i.e. vector
    static const bool IsScalar = false;
    static const bool IsVector = true;
    static const int  Size = T::SIZE;
    using ElementType = typename T::ValueType;
    static ElementType scalar(const T& v) { return v.length(); }
};

// ----------------------------> FloatTraits <--------------------------------------

template<typename T, int = sizeof(typename TensorTraits<T>::ElementType)>
struct FloatTraits
{
    using FloatType = float;
};

template<typename T>
struct FloatTraits<T, 8>
{
    using FloatType = double;
};

template<>
struct FloatTraits<bool, 1>
{
    using FloatType = bool;
};

template<>
struct FloatTraits<ValueIndex, 1> // size of empty class in C++ is 1 byte and not 0 byte
{
    using FloatType = uint64_t;
};

template<>
struct FloatTraits<ValueIndexMask, 1> // size of empty class in C++ is 1 byte and not 0 byte
{
    using FloatType = uint64_t;
};

template<>
struct FloatTraits<ValueOnIndex, 1> // size of empty class in C++ is 1 byte and not 0 byte
{
    using FloatType = uint64_t;
};

template<>
struct FloatTraits<ValueOnIndexMask, 1> // size of empty class in C++ is 1 byte and not 0 byte
{
    using FloatType = uint64_t;
};

template<>
struct FloatTraits<ValueMask, 1> // size of empty class in C++ is 1 byte and not 0 byte
{
    using FloatType = bool;
};

template<>
struct FloatTraits<Point, 1> // size of empty class in C++ is 1 byte and not 0 byte
{
    using FloatType = double;
};

// ----------------------------> mapping BuildType -> GridType <--------------------------------------

/// @brief Maps from a templated build type to a GridType enum
template<typename BuildT>
__hostdev__ inline GridType mapToGridType()
{
    if constexpr(is_same<BuildT, float>::value) { // resolved at compile-time
        return GridType::Float;
    } else if constexpr(is_same<BuildT, double>::value) {
        return GridType::Double;
    } else if constexpr(is_same<BuildT, int16_t>::value) {
        return GridType::Int16;
    } else if constexpr(is_same<BuildT, int32_t>::value) {
        return GridType::Int32;
    } else if constexpr(is_same<BuildT, int64_t>::value) {
        return GridType::Int64;
    } else if constexpr(is_same<BuildT, Vec3f>::value) {
        return GridType::Vec3f;
    } else if constexpr(is_same<BuildT, Vec3d>::value) {
        return GridType::Vec3d;
    } else if constexpr(is_same<BuildT, uint32_t>::value) {
        return GridType::UInt32;
    } else if constexpr(is_same<BuildT, ValueMask>::value) {
        return GridType::Mask;
    } else if constexpr(is_same<BuildT, Half>::value) {
        return GridType::Half;
    } else if constexpr(is_same<BuildT, ValueIndex>::value) {
        return GridType::Index;
    } else if constexpr(is_same<BuildT, ValueOnIndex>::value) {
        return GridType::OnIndex;
    } else if constexpr(is_same<BuildT, ValueIndexMask>::value) {
        return GridType::IndexMask;
    } else if constexpr(is_same<BuildT, ValueOnIndexMask>::value) {
        return GridType::OnIndexMask;
    } else if constexpr(is_same<BuildT, bool>::value) {
        return GridType::Boolean;
    } else if constexpr(is_same<BuildT, Rgba8>::value) {
        return GridType::RGBA8;
    } else if (is_same<BuildT, Fp4>::value) {
        return GridType::Fp4;
    } else if constexpr(is_same<BuildT, Fp8>::value) {
        return GridType::Fp8;
    } else if constexpr(is_same<BuildT, Fp16>::value) {
        return GridType::Fp16;
    } else if constexpr(is_same<BuildT, FpN>::value) {
        return GridType::FpN;
    } else if constexpr(is_same<BuildT, Vec4f>::value) {
        return GridType::Vec4f;
    } else if constexpr(is_same<BuildT, Vec4d>::value) {
        return GridType::Vec4d;
    } else if (is_same<BuildT, Point>::value) {
        return GridType::PointIndex;
    } else if constexpr(is_same<BuildT, Vec3u8>::value) {
        return GridType::Vec3u8;
    } else if constexpr(is_same<BuildT, Vec3u16>::value) {
        return GridType::Vec3u16;
    }
    return GridType::Unknown;
}

// ----------------------------> mapping BuildType -> GridClass <--------------------------------------

/// @brief Maps from a templated build type to a GridClass enum
template<typename BuildT>
__hostdev__ inline GridClass mapToGridClass(GridClass defaultClass = GridClass::Unknown)
{
    if (is_same<BuildT, ValueMask>::value) {
        return GridClass::Topology;
    } else if (BuildTraits<BuildT>::is_index) {
        return GridClass::IndexGrid;
    } else if (is_same<BuildT, Rgba8>::value) {
        return GridClass::VoxelVolume;
    } else if (is_same<BuildT, Point>::value) {
        return GridClass::PointIndex;
    }
    return defaultClass;
}

// ----------------------------> matMult <--------------------------------------

/// @brief Multiply a 3x3 matrix and a 3d vector using 32bit floating point arithmetics
/// @note This corresponds to a linear mapping, e.g. scaling, rotation etc.
/// @tparam Vec3T Template type of the input and output 3d vectors
/// @param mat pointer to an array of floats with the 3x3 matrix
/// @param xyz input vector to be multiplied by the matrix
/// @return result of matrix-vector multiplication, i.e. mat x xyz
template<typename Vec3T>
__hostdev__ inline Vec3T matMult(const float* mat, const Vec3T& xyz)
{
    return Vec3T(fmaf(static_cast<float>(xyz[0]), mat[0], fmaf(static_cast<float>(xyz[1]), mat[1], static_cast<float>(xyz[2]) * mat[2])),
                 fmaf(static_cast<float>(xyz[0]), mat[3], fmaf(static_cast<float>(xyz[1]), mat[4], static_cast<float>(xyz[2]) * mat[5])),
                 fmaf(static_cast<float>(xyz[0]), mat[6], fmaf(static_cast<float>(xyz[1]), mat[7], static_cast<float>(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops
}

/// @brief Multiply a 3x3 matrix and a 3d vector using 64bit floating point arithmetics
/// @note This corresponds to a linear mapping, e.g. scaling, rotation etc.
/// @tparam Vec3T Template type of the input and output 3d vectors
/// @param mat pointer to an array of floats with the 3x3 matrix
/// @param xyz input vector to be multiplied by the matrix
/// @return result of matrix-vector multiplication, i.e. mat x xyz
template<typename Vec3T>
__hostdev__ inline Vec3T matMult(const double* mat, const Vec3T& xyz)
{
    return Vec3T(fma(static_cast<double>(xyz[0]), mat[0], fma(static_cast<double>(xyz[1]), mat[1], static_cast<double>(xyz[2]) * mat[2])),
                 fma(static_cast<double>(xyz[0]), mat[3], fma(static_cast<double>(xyz[1]), mat[4], static_cast<double>(xyz[2]) * mat[5])),
                 fma(static_cast<double>(xyz[0]), mat[6], fma(static_cast<double>(xyz[1]), mat[7], static_cast<double>(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops
}

/// @brief Multiply a 3x3 matrix to a 3d vector and add another 3d vector using 32bit floating point arithmetics
/// @note This corresponds to an affine transformation, i.e a linear mapping followed by a translation. e.g. scale/rotation and translation
/// @tparam Vec3T Template type of the input and output 3d vectors
/// @param mat pointer to an array of floats with the 3x3 matrix
/// @param vec 3d vector to be added AFTER the matrix multiplication
/// @param xyz input vector to be multiplied by the matrix and a translated by @c vec
/// @return result of affine transformation, i.e. (mat x xyz) + vec
template<typename Vec3T>
__hostdev__ inline Vec3T matMult(const float* mat, const float* vec, const Vec3T& xyz)
{
    return Vec3T(fmaf(static_cast<float>(xyz[0]), mat[0], fmaf(static_cast<float>(xyz[1]), mat[1], fmaf(static_cast<float>(xyz[2]), mat[2], vec[0]))),
                 fmaf(static_cast<float>(xyz[0]), mat[3], fmaf(static_cast<float>(xyz[1]), mat[4], fmaf(static_cast<float>(xyz[2]), mat[5], vec[1]))),
                 fmaf(static_cast<float>(xyz[0]), mat[6], fmaf(static_cast<float>(xyz[1]), mat[7], fmaf(static_cast<float>(xyz[2]), mat[8], vec[2])))); // 9 fmaf = 9 flops
}

/// @brief Multiply a 3x3 matrix to a 3d vector and add another 3d vector using 64bit floating point arithmetics
/// @note This corresponds to an affine transformation, i.e a linear mapping followed by a translation. e.g. scale/rotation and translation
/// @tparam Vec3T Template type of the input and output 3d vectors
/// @param mat pointer to an array of floats with the 3x3 matrix
/// @param vec 3d vector to be added AFTER the matrix multiplication
/// @param xyz input vector to be multiplied by the matrix and a translated by @c vec
/// @return result of affine transformation, i.e. (mat x xyz) + vec
template<typename Vec3T>
__hostdev__ inline Vec3T matMult(const double* mat, const double* vec, const Vec3T& xyz)
{
    return Vec3T(fma(static_cast<double>(xyz[0]), mat[0], fma(static_cast<double>(xyz[1]), mat[1], fma(static_cast<double>(xyz[2]), mat[2], vec[0]))),
                 fma(static_cast<double>(xyz[0]), mat[3], fma(static_cast<double>(xyz[1]), mat[4], fma(static_cast<double>(xyz[2]), mat[5], vec[1]))),
                 fma(static_cast<double>(xyz[0]), mat[6], fma(static_cast<double>(xyz[1]), mat[7], fma(static_cast<double>(xyz[2]), mat[8], vec[2])))); // 9 fma = 9 flops
}

/// @brief Multiply the transposed of a 3x3 matrix and a 3d vector using 32bit floating point arithmetics
/// @note This corresponds to an inverse linear mapping, e.g. inverse scaling, inverse rotation etc.
/// @tparam Vec3T Template type of the input and output 3d vectors
/// @param mat pointer to an array of floats with the 3x3 matrix
/// @param xyz input vector to be multiplied by the transposed matrix
/// @return result of matrix-vector multiplication, i.e. mat^T x xyz
template<typename Vec3T>
__hostdev__ inline Vec3T matMultT(const float* mat, const Vec3T& xyz)
{
    return Vec3T(fmaf(static_cast<float>(xyz[0]), mat[0], fmaf(static_cast<float>(xyz[1]), mat[3], static_cast<float>(xyz[2]) * mat[6])),
                 fmaf(static_cast<float>(xyz[0]), mat[1], fmaf(static_cast<float>(xyz[1]), mat[4], static_cast<float>(xyz[2]) * mat[7])),
                 fmaf(static_cast<float>(xyz[0]), mat[2], fmaf(static_cast<float>(xyz[1]), mat[5], static_cast<float>(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops
}

/// @brief Multiply the transposed of a 3x3 matrix and a 3d vector using 64bit floating point arithmetics
/// @note This corresponds to an inverse linear mapping, e.g. inverse scaling, inverse rotation etc.
/// @tparam Vec3T Template type of the input and output 3d vectors
/// @param mat pointer to an array of floats with the 3x3 matrix
/// @param xyz input vector to be multiplied by the transposed matrix
/// @return result of matrix-vector multiplication, i.e. mat^T x xyz
template<typename Vec3T>
__hostdev__ inline Vec3T matMultT(const double* mat, const Vec3T& xyz)
{
    return Vec3T(fma(static_cast<double>(xyz[0]), mat[0], fma(static_cast<double>(xyz[1]), mat[3], static_cast<double>(xyz[2]) * mat[6])),
                 fma(static_cast<double>(xyz[0]), mat[1], fma(static_cast<double>(xyz[1]), mat[4], static_cast<double>(xyz[2]) * mat[7])),
                 fma(static_cast<double>(xyz[0]), mat[2], fma(static_cast<double>(xyz[1]), mat[5], static_cast<double>(xyz[2]) * mat[8]))); // 6 fmaf + 3 mult = 9 flops
}

template<typename Vec3T>
__hostdev__ inline Vec3T matMultT(const float* mat, const float* vec, const Vec3T& xyz)
{
    return Vec3T(fmaf(static_cast<float>(xyz[0]), mat[0], fmaf(static_cast<float>(xyz[1]), mat[3], fmaf(static_cast<float>(xyz[2]), mat[6], vec[0]))),
                 fmaf(static_cast<float>(xyz[0]), mat[1], fmaf(static_cast<float>(xyz[1]), mat[4], fmaf(static_cast<float>(xyz[2]), mat[7], vec[1]))),
                 fmaf(static_cast<float>(xyz[0]), mat[2], fmaf(static_cast<float>(xyz[1]), mat[5], fmaf(static_cast<float>(xyz[2]), mat[8], vec[2])))); // 9 fmaf = 9 flops
}

template<typename Vec3T>
__hostdev__ inline Vec3T matMultT(const double* mat, const double* vec, const Vec3T& xyz)
{
    return Vec3T(fma(static_cast<double>(xyz[0]), mat[0], fma(static_cast<double>(xyz[1]), mat[3], fma(static_cast<double>(xyz[2]), mat[6], vec[0]))),
                 fma(static_cast<double>(xyz[0]), mat[1], fma(static_cast<double>(xyz[1]), mat[4], fma(static_cast<double>(xyz[2]), mat[7], vec[1]))),
                 fma(static_cast<double>(xyz[0]), mat[2], fma(static_cast<double>(xyz[1]), mat[5], fma(static_cast<double>(xyz[2]), mat[8], vec[2])))); // 9 fma = 9 flops
}

// ----------------------------> BBox <-------------------------------------

// Base-class for static polymorphism (cannot be constructed directly)
template<typename Vec3T>
struct BaseBBox
{
    Vec3T                    mCoord[2];
    __hostdev__ bool         operator==(const BaseBBox& rhs) const { return mCoord[0] == rhs.mCoord[0] && mCoord[1] == rhs.mCoord[1]; };
    __hostdev__ bool         operator!=(const BaseBBox& rhs) const { return mCoord[0] != rhs.mCoord[0] || mCoord[1] != rhs.mCoord[1]; };
    __hostdev__ const Vec3T& operator[](int i) const { return mCoord[i]; }
    __hostdev__ Vec3T&       operator[](int i) { return mCoord[i]; }
    __hostdev__ Vec3T&       min() { return mCoord[0]; }
    __hostdev__ Vec3T&       max() { return mCoord[1]; }
    __hostdev__ const Vec3T& min() const { return mCoord[0]; }
    __hostdev__ const Vec3T& max() const { return mCoord[1]; }
    __hostdev__ BaseBBox&    translate(const Vec3T& xyz)
    {
        mCoord[0] += xyz;
        mCoord[1] += xyz;
        return *this;
    }
    /// @brief Expand this bounding box to enclose point @c xyz.
    __hostdev__ BaseBBox& expand(const Vec3T& xyz)
    {
        mCoord[0].minComponent(xyz);
        mCoord[1].maxComponent(xyz);
        return *this;
    }

    /// @brief Expand this bounding box to enclose the given bounding box.
    __hostdev__ BaseBBox& expand(const BaseBBox& bbox)
    {
        mCoord[0].minComponent(bbox[0]);
        mCoord[1].maxComponent(bbox[1]);
        return *this;
    }

    /// @brief Intersect this bounding box with the given bounding box.
    __hostdev__ BaseBBox& intersect(const BaseBBox& bbox)
    {
        mCoord[0].maxComponent(bbox[0]);
        mCoord[1].minComponent(bbox[1]);
        return *this;
    }

    //__hostdev__ BaseBBox expandBy(typename Vec3T::ValueType padding) const
    //{
    //    return BaseBBox(mCoord[0].offsetBy(-padding),mCoord[1].offsetBy(padding));
    //}
    __hostdev__ bool isInside(const Vec3T& xyz)
    {
        if (xyz[0] < mCoord[0][0] || xyz[1] < mCoord[0][1] || xyz[2] < mCoord[0][2])
            return false;
        if (xyz[0] > mCoord[1][0] || xyz[1] > mCoord[1][1] || xyz[2] > mCoord[1][2])
            return false;
        return true;
    }

protected:
    __hostdev__ BaseBBox() {}
    __hostdev__ BaseBBox(const Vec3T& min, const Vec3T& max)
        : mCoord{min, max}
    {
    }
}; // BaseBBox

template<typename Vec3T, bool = is_floating_point<typename Vec3T::ValueType>::value>
struct BBox;

/// @brief Partial template specialization for floating point coordinate types.
///
/// @note Min is inclusive and max is exclusive. If min = max the dimension of
///       the bounding box is zero and therefore it is also empty.
template<typename Vec3T>
struct BBox<Vec3T, true> : public BaseBBox<Vec3T>
{
    using Vec3Type = Vec3T;
    using ValueType = typename Vec3T::ValueType;
    static_assert(is_floating_point<ValueType>::value, "Expected a floating point coordinate type");
    using BaseT = BaseBBox<Vec3T>;
    using BaseT::mCoord;
    /// @brief Default construction sets BBox to an empty bbox
    __hostdev__ BBox()
        : BaseT(Vec3T( Maximum<typename Vec3T::ValueType>::value()),
                Vec3T(-Maximum<typename Vec3T::ValueType>::value()))
    {
    }
    __hostdev__ BBox(const Vec3T& min, const Vec3T& max)
        : BaseT(min, max)
    {
    }
    __hostdev__ BBox(const Coord& min, const Coord& max)
        : BaseT(Vec3T(ValueType(min[0]), ValueType(min[1]), ValueType(min[2])),
                Vec3T(ValueType(max[0] + 1), ValueType(max[1] + 1), ValueType(max[2] + 1)))
    {
    }
    __hostdev__ static BBox createCube(const Coord& min, typename Coord::ValueType dim)
    {
        return BBox(min, min.offsetBy(dim));
    }

    __hostdev__ BBox(const BaseBBox<Coord>& bbox)
        : BBox(bbox[0], bbox[1])
    {
    }
    __hostdev__ bool  empty() const { return mCoord[0][0] >= mCoord[1][0] ||
                                             mCoord[0][1] >= mCoord[1][1] ||
                                             mCoord[0][2] >= mCoord[1][2]; }
    __hostdev__ operator bool() const { return mCoord[0][0] < mCoord[1][0] &&
                                               mCoord[0][1] < mCoord[1][1] &&
                                               mCoord[0][2] < mCoord[1][2]; }
    __hostdev__ Vec3T dim() const { return *this ? this->max() - this->min() : Vec3T(0); }
    __hostdev__ bool  isInside(const Vec3T& p) const
    {
        return p[0] > mCoord[0][0] && p[1] > mCoord[0][1] && p[2] > mCoord[0][2] &&
               p[0] < mCoord[1][0] && p[1] < mCoord[1][1] && p[2] < mCoord[1][2];
    }

}; // BBox<Vec3T, true>

/// @brief Partial template specialization for integer coordinate types
///
/// @note Both min and max are INCLUDED in the bbox so dim = max - min + 1. So,
///       if min = max the bounding box contains exactly one point and dim = 1!
template<typename CoordT>
struct BBox<CoordT, false> : public BaseBBox<CoordT>
{
    static_assert(is_same<int, typename CoordT::ValueType>::value, "Expected \"int\" coordinate type");
    using BaseT = BaseBBox<CoordT>;
    using BaseT::mCoord;
    /// @brief Iterator over the domain covered by a BBox
    /// @details z is the fastest-moving coordinate.
    class Iterator
    {
        const BBox& mBBox;
        CoordT      mPos;

    public:
        __hostdev__ Iterator(const BBox& b)
            : mBBox(b)
            , mPos(b.min())
        {
        }
        __hostdev__ Iterator(const BBox& b, const Coord& p)
            : mBBox(b)
            , mPos(p)
        {
        }
        __hostdev__ Iterator& operator++()
        {
            if (mPos[2] < mBBox[1][2]) { // this is the most common case
                ++mPos[2];// increment z
            } else if (mPos[1] < mBBox[1][1]) {
                mPos[2] = mBBox[0][2];// reset z
                ++mPos[1];// increment y
            } else if (mPos[0] <= mBBox[1][0]) {
                mPos[2] = mBBox[0][2];// reset z
                mPos[1] = mBBox[0][1];// reset y
                ++mPos[0];// increment x
            }
            return *this;
        }
        __hostdev__ Iterator operator++(int)
        {
            auto tmp = *this;
            ++(*this);
            return tmp;
        }
        __hostdev__ bool operator==(const Iterator& rhs) const
        {
            NANOVDB_ASSERT(mBBox == rhs.mBBox);
            return mPos == rhs.mPos;
        }
        __hostdev__ bool operator!=(const Iterator& rhs) const
        {
            NANOVDB_ASSERT(mBBox == rhs.mBBox);
            return mPos != rhs.mPos;
        }
        __hostdev__ bool operator<(const Iterator& rhs) const
        {
            NANOVDB_ASSERT(mBBox == rhs.mBBox);
            return mPos < rhs.mPos;
        }
        __hostdev__ bool operator<=(const Iterator& rhs) const
        {
            NANOVDB_ASSERT(mBBox == rhs.mBBox);
            return mPos <= rhs.mPos;
        }
        /// @brief Return @c true if the iterator still points to a valid coordinate.
        __hostdev__ operator bool() const { return mPos <= mBBox[1]; }
        __hostdev__ const CoordT& operator*() const { return mPos; }
    }; // Iterator
    __hostdev__ Iterator begin() const { return Iterator{*this}; }
    __hostdev__ Iterator end()   const { return Iterator{*this, CoordT(mCoord[1][0]+1, mCoord[0][1], mCoord[0][2])}; }
    __hostdev__          BBox()
        : BaseT(CoordT::max(), CoordT::min())
    {
    }
    __hostdev__ BBox(const CoordT& min, const CoordT& max)
        : BaseT(min, max)
    {
    }

    template<typename SplitT>
    __hostdev__ BBox(BBox& other, const SplitT&)
        : BaseT(other.mCoord[0], other.mCoord[1])
    {
        NANOVDB_ASSERT(this->is_divisible());
        const int n = MaxIndex(this->dim());
        mCoord[1][n] = (mCoord[0][n] + mCoord[1][n]) >> 1;
        other.mCoord[0][n] = mCoord[1][n] + 1;
    }

    __hostdev__ static BBox createCube(const CoordT& min, typename CoordT::ValueType dim)
    {
        return BBox(min, min.offsetBy(dim - 1));
    }

    __hostdev__ static BBox createCube(typename CoordT::ValueType min, typename CoordT::ValueType max)
    {
        return BBox(CoordT(min), CoordT(max));
    }

    __hostdev__ bool is_divisible() const { return mCoord[0][0] < mCoord[1][0] &&
                                                   mCoord[0][1] < mCoord[1][1] &&
                                                   mCoord[0][2] < mCoord[1][2]; }
    /// @brief Return true if this bounding box is empty, e.g. uninitialized
    __hostdev__ bool     empty() const { return mCoord[0][0] > mCoord[1][0] ||
                                                mCoord[0][1] > mCoord[1][1] ||
                                                mCoord[0][2] > mCoord[1][2]; }
    /// @brief Convert this BBox to boolean true if it is not empty
    __hostdev__ operator bool() const { return mCoord[0][0] <= mCoord[1][0] &&
                                               mCoord[0][1] <= mCoord[1][1] &&
                                               mCoord[0][2] <= mCoord[1][2]; }
    __hostdev__ CoordT   dim() const { return *this ? this->max() - this->min() + Coord(1) : Coord(0); }
    __hostdev__ uint64_t volume() const
    {
        auto d = this->dim();
        return uint64_t(d[0]) * uint64_t(d[1]) * uint64_t(d[2]);
    }
    __hostdev__ bool isInside(const CoordT& p) const { return !(CoordT::lessThan(p, this->min()) || CoordT::lessThan(this->max(), p)); }
    /// @brief Return @c true if the given bounding box is inside this bounding box.
    __hostdev__ bool isInside(const BBox& b) const
    {
        return !(CoordT::lessThan(b.min(), this->min()) || CoordT::lessThan(this->max(), b.max()));
    }

    /// @brief Return @c true if the given bounding box overlaps with this bounding box.
    __hostdev__ bool hasOverlap(const BBox& b) const
    {
        return !(CoordT::lessThan(this->max(), b.min()) || CoordT::lessThan(b.max(), this->min()));
    }

    /// @warning This converts a CoordBBox into a floating-point bounding box which implies that max += 1 !
    template<typename RealT = double>
    __hostdev__ BBox<Vec3<RealT>> asReal() const
    {
        static_assert(is_floating_point<RealT>::value, "CoordBBox::asReal: Expected a floating point coordinate");
        return BBox<Vec3<RealT>>(Vec3<RealT>(RealT(mCoord[0][0]), RealT(mCoord[0][1]), RealT(mCoord[0][2])),
                                 Vec3<RealT>(RealT(mCoord[1][0] + 1), RealT(mCoord[1][1] + 1), RealT(mCoord[1][2] + 1)));
    }
    /// @brief Return a new instance that is expanded by the specified padding.
    __hostdev__ BBox expandBy(typename CoordT::ValueType padding) const
    {
        return BBox(mCoord[0].offsetBy(-padding), mCoord[1].offsetBy(padding));
    }

    /// @brief  @brief transform this coordinate bounding box by the specified map
    /// @param map mapping of index to world coordinates
    /// @return world bounding box
    template<typename Map>
    __hostdev__ BBox<Vec3d> transform(const Map& map) const
    {
        const Vec3d tmp = map.applyMap(Vec3d(mCoord[0][0], mCoord[0][1], mCoord[0][2]));
        BBox<Vec3d> bbox(tmp, tmp);
        bbox.expand(map.applyMap(Vec3d(mCoord[0][0], mCoord[0][1], mCoord[1][2])));
        bbox.expand(map.applyMap(Vec3d(mCoord[0][0], mCoord[1][1], mCoord[0][2])));
        bbox.expand(map.applyMap(Vec3d(mCoord[1][0], mCoord[0][1], mCoord[0][2])));
        bbox.expand(map.applyMap(Vec3d(mCoord[1][0], mCoord[1][1], mCoord[0][2])));
        bbox.expand(map.applyMap(Vec3d(mCoord[1][0], mCoord[0][1], mCoord[1][2])));
        bbox.expand(map.applyMap(Vec3d(mCoord[0][0], mCoord[1][1], mCoord[1][2])));
        bbox.expand(map.applyMap(Vec3d(mCoord[1][0], mCoord[1][1], mCoord[1][2])));
        return bbox;
    }

#if defined(__CUDACC__) // the following functions only run on the GPU!
    __device__ inline BBox& expandAtomic(const CoordT& ijk)
    {
        mCoord[0].minComponentAtomic(ijk);
        mCoord[1].maxComponentAtomic(ijk);
        return *this;
    }
    __device__ inline BBox& expandAtomic(const BBox& bbox)
    {
        mCoord[0].minComponentAtomic(bbox[0]);
        mCoord[1].maxComponentAtomic(bbox[1]);
        return *this;
    }
    __device__ inline BBox& intersectAtomic(const BBox& bbox)
    {
        mCoord[0].maxComponentAtomic(bbox[0]);
        mCoord[1].minComponentAtomic(bbox[1]);
        return *this;
    }
#endif
}; // BBox<CoordT, false>

using CoordBBox = BBox<Coord>;
using BBoxR = BBox<Vec3d>;

// -------------------> Find lowest and highest bit in a word <----------------------------

/// @brief Returns the index of the lowest, i.e. least significant, on bit in the specified 32 bit word
///
/// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
NANOVDB_HOSTDEV_DISABLE_WARNING
__hostdev__ static inline uint32_t FindLowestOn(uint32_t v)
{
    NANOVDB_ASSERT(v);
#if (defined(__CUDA_ARCH__) || defined(__HIP__)) && defined(NANOVDB_USE_INTRINSICS)
    return __ffs(v) - 1; // one based indexing
#elif defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
    unsigned long index;
    _BitScanForward(&index, v);
    return static_cast<uint32_t>(index);
#elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
    return static_cast<uint32_t>(__builtin_ctzl(v));
#else
    //NANO_WARNING("Using software implementation for FindLowestOn(uint32_t v)")
    static const unsigned char DeBruijn[32] = {
        0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9};
// disable unary minus on unsigned warning
#if defined(_MSC_VER) && !defined(__NVCC__)
#pragma warning(push)
#pragma warning(disable : 4146)
#endif
    return DeBruijn[uint32_t((v & -v) * 0x077CB531U) >> 27];
#if defined(_MSC_VER) && !defined(__NVCC__)
#pragma warning(pop)
#endif

#endif
}

/// @brief Returns the index of the highest, i.e. most significant, on bit in the specified 32 bit word
///
/// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
NANOVDB_HOSTDEV_DISABLE_WARNING
__hostdev__ static inline uint32_t FindHighestOn(uint32_t v)
{
    NANOVDB_ASSERT(v);
#if (defined(__CUDA_ARCH__) || defined(__HIP__)) && defined(NANOVDB_USE_INTRINSICS)
    return sizeof(uint32_t) * 8 - 1 - __clz(v); // Return the number of consecutive high-order zero bits in a 32-bit integer.
#elif defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
    unsigned long index;
    _BitScanReverse(&index, v);
    return static_cast<uint32_t>(index);
#elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
    return sizeof(unsigned long) * 8 - 1 - __builtin_clzl(v);
#else
    //NANO_WARNING("Using software implementation for FindHighestOn(uint32_t)")
    static const unsigned char DeBruijn[32] = {
        0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
        8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31};
    v |= v >> 1; // first round down to one less than a power of 2
    v |= v >> 2;
    v |= v >> 4;
    v |= v >> 8;
    v |= v >> 16;
    return DeBruijn[uint32_t(v * 0x07C4ACDDU) >> 27];
#endif
}

/// @brief Returns the index of the lowest, i.e. least significant, on bit in the specified 64 bit word
///
/// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
NANOVDB_HOSTDEV_DISABLE_WARNING
__hostdev__ static inline uint32_t FindLowestOn(uint64_t v)
{
    NANOVDB_ASSERT(v);
#if (defined(__CUDA_ARCH__) || defined(__HIP__)) && defined(NANOVDB_USE_INTRINSICS)
    return __ffsll(static_cast<unsigned long long int>(v)) - 1; // one based indexing
#elif defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
    unsigned long index;
    _BitScanForward64(&index, v);
    return static_cast<uint32_t>(index);
#elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
    return static_cast<uint32_t>(__builtin_ctzll(v));
#else
    //NANO_WARNING("Using software implementation for FindLowestOn(uint64_t)")
    static const unsigned char DeBruijn[64] = {
        0,   1,  2, 53,  3,  7, 54, 27, 4,  38, 41,  8, 34, 55, 48, 28,
        62,  5, 39, 46, 44, 42, 22,  9, 24, 35, 59, 56, 49, 18, 29, 11,
        63, 52,  6, 26, 37, 40, 33, 47, 61, 45, 43, 21, 23, 58, 17, 10,
        51, 25, 36, 32, 60, 20, 57, 16, 50, 31, 19, 15, 30, 14, 13, 12,
    };
// disable unary minus on unsigned warning
#if defined(_MSC_VER) && !defined(__NVCC__)
#pragma warning(push)
#pragma warning(disable : 4146)
#endif
    return DeBruijn[uint64_t((v & -v) * UINT64_C(0x022FDD63CC95386D)) >> 58];
#if defined(_MSC_VER) && !defined(__NVCC__)
#pragma warning(pop)
#endif

#endif
}

/// @brief Returns the index of the highest, i.e. most significant, on bit in the specified 64 bit word
///
/// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
NANOVDB_HOSTDEV_DISABLE_WARNING
__hostdev__ static inline uint32_t FindHighestOn(uint64_t v)
{
    NANOVDB_ASSERT(v);
#if (defined(__CUDA_ARCH__) || defined(__HIP__)) && defined(NANOVDB_USE_INTRINSICS)
    return sizeof(unsigned long) * 8 - 1 - __clzll(static_cast<unsigned long long int>(v));
#elif defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
    unsigned long index;
    _BitScanReverse64(&index, v);
    return static_cast<uint32_t>(index);
#elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
    return sizeof(unsigned long) * 8 - 1 - __builtin_clzll(v);
#else
    const uint32_t* p = reinterpret_cast<const uint32_t*>(&v);
    return p[1] ? 32u + FindHighestOn(p[1]) : FindHighestOn(p[0]);
#endif
}

// ----------------------------> CountOn <--------------------------------------

/// @return Number of bits that are on in the specified 64-bit word
NANOVDB_HOSTDEV_DISABLE_WARNING
__hostdev__ inline uint32_t CountOn(uint64_t v)
{
#if (defined(__CUDA_ARCH__) || defined(__HIP__)) && defined(NANOVDB_USE_INTRINSICS)
    //#warning Using popcll for CountOn
    return __popcll(v);
// __popcnt64 intrinsic support was added in VS 2019 16.8
#elif defined(_MSC_VER) && defined(_M_X64) && (_MSC_VER >= 1928) && defined(NANOVDB_USE_INTRINSICS)
    //#warning Using popcnt64 for CountOn
    return uint32_t(__popcnt64(v));
#elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
    //#warning Using builtin_popcountll for CountOn
    return __builtin_popcountll(v);
#else // use software implementation
    //NANO_WARNING("Using software implementation for CountOn")
    v = v - ((v >> 1) & uint64_t(0x5555555555555555));
    v = (v & uint64_t(0x3333333333333333)) + ((v >> 2) & uint64_t(0x3333333333333333));
    return (((v + (v >> 4)) & uint64_t(0xF0F0F0F0F0F0F0F)) * uint64_t(0x101010101010101)) >> 56;
#endif
}

//  ----------------------------> BitFlags <--------------------------------------

template<int N>
struct BitArray;
template<>
struct BitArray<8>
{
    uint8_t mFlags{0};
};
template<>
struct BitArray<16>
{
    uint16_t mFlags{0};
};
template<>
struct BitArray<32>
{
    uint32_t mFlags{0};
};
template<>
struct BitArray<64>
{
    uint64_t mFlags{0};
};

template<int N>
class BitFlags : public BitArray<N>
{
protected:
    using BitArray<N>::mFlags;

public:
    using Type = decltype(mFlags);
    BitFlags() {}
    BitFlags(std::initializer_list<uint8_t> list)
    {
        for (auto bit : list)
            mFlags |= static_cast<Type>(1 << bit);
    }
    template<typename MaskT>
    BitFlags(std::initializer_list<MaskT> list)
    {
        for (auto mask : list)
            mFlags |= static_cast<Type>(mask);
    }
    __hostdev__ Type  data() const { return mFlags; }
    __hostdev__ Type& data() { return mFlags; }
    __hostdev__ void  initBit(std::initializer_list<uint8_t> list)
    {
        mFlags = 0u;
        for (auto bit : list)
            mFlags |= static_cast<Type>(1 << bit);
    }
    template<typename MaskT>
    __hostdev__ void initMask(std::initializer_list<MaskT> list)
    {
        mFlags = 0u;
        for (auto mask : list)
            mFlags |= static_cast<Type>(mask);
    }
    //__hostdev__ Type& data() { return mFlags; }
    //__hostdev__ Type data() const { return mFlags; }
    __hostdev__ Type getFlags() const { return mFlags & (static_cast<Type>(GridFlags::End) - 1u); } // mask out everything except relevant bits

    __hostdev__ void setOn() { mFlags = ~Type(0u); }
    __hostdev__ void setOff() { mFlags = Type(0u); }

    __hostdev__ void setBitOn(uint8_t bit) { mFlags |= static_cast<Type>(1 << bit); }
    __hostdev__ void setBitOff(uint8_t bit) { mFlags &= ~static_cast<Type>(1 << bit); }

    __hostdev__ void setBitOn(std::initializer_list<uint8_t> list)
    {
        for (auto bit : list)
            mFlags |= static_cast<Type>(1 << bit);
    }
    __hostdev__ void setBitOff(std::initializer_list<uint8_t> list)
    {
        for (auto bit : list)
            mFlags &= ~static_cast<Type>(1 << bit);
    }

    template<typename MaskT>
    __hostdev__ void setMaskOn(MaskT mask) { mFlags |= static_cast<Type>(mask); }
    template<typename MaskT>
    __hostdev__ void setMaskOff(MaskT mask) { mFlags &= ~static_cast<Type>(mask); }

    template<typename MaskT>
    __hostdev__ void setMaskOn(std::initializer_list<MaskT> list)
    {
        for (auto mask : list)
            mFlags |= static_cast<Type>(mask);
    }
    template<typename MaskT>
    __hostdev__ void setMaskOff(std::initializer_list<MaskT> list)
    {
        for (auto mask : list)
            mFlags &= ~static_cast<Type>(mask);
    }

    __hostdev__ void setBit(uint8_t bit, bool on) { on ? this->setBitOn(bit) : this->setBitOff(bit); }
    template<typename MaskT>
    __hostdev__ void setMask(MaskT mask, bool on) { on ? this->setMaskOn(mask) : this->setMaskOff(mask); }

    __hostdev__ bool isOn() const { return mFlags == ~Type(0u); }
    __hostdev__ bool isOff() const { return mFlags == Type(0u); }
    __hostdev__ bool isBitOn(uint8_t bit) const { return 0 != (mFlags & static_cast<Type>(1 << bit)); }
    __hostdev__ bool isBitOff(uint8_t bit) const { return 0 == (mFlags & static_cast<Type>(1 << bit)); }
    template<typename MaskT>
    __hostdev__ bool isMaskOn(MaskT mask) const { return 0 != (mFlags & static_cast<Type>(mask)); }
    template<typename MaskT>
    __hostdev__ bool isMaskOff(MaskT mask) const { return 0 == (mFlags & static_cast<Type>(mask)); }
    /// @brief return true if any of the masks in the list are on
    template<typename MaskT>
    __hostdev__ bool isMaskOn(std::initializer_list<MaskT> list) const
    {
        for (auto mask : list)
            if (0 != (mFlags & static_cast<Type>(mask)))
                return true;
        return false;
    }
    /// @brief return true if any of the masks in the list are off
    template<typename MaskT>
    __hostdev__ bool isMaskOff(std::initializer_list<MaskT> list) const
    {
        for (auto mask : list)
            if (0 == (mFlags & static_cast<Type>(mask)))
                return true;
        return false;
    }
    /// @brief required for backwards compatibility
    __hostdev__ BitFlags& operator=(Type n)
    {
        mFlags = n;
        return *this;
    }
}; // BitFlags<N>

// ----------------------------> Mask <--------------------------------------

/// @brief Bit-mask to encode active states and facilitate sequential iterators
/// and a fast codec for I/O compression.
template<uint32_t LOG2DIM>
class Mask
{
public:
    static constexpr uint32_t SIZE = 1U << (3 * LOG2DIM); // Number of bits in mask
    static constexpr uint32_t WORD_COUNT = SIZE >> 6; // Number of 64 bit words

    /// @brief Return the memory footprint in bytes of this Mask
    __hostdev__ static size_t memUsage() { return sizeof(Mask); }

    /// @brief Return the number of bits available in this Mask
    __hostdev__ static uint32_t bitCount() { return SIZE; }

    /// @brief Return the number of machine words used by this Mask
    __hostdev__ static uint32_t wordCount() { return WORD_COUNT; }

    /// @brief Return the total number of set bits in this Mask
    __hostdev__ uint32_t countOn() const
    {
        uint32_t sum = 0;
        for (const uint64_t *w = mWords, *q = w + WORD_COUNT; w != q; ++w)
            sum += CountOn(*w);
        return sum;
    }

    /// @brief Return the number of lower set bits in mask up to but excluding the i'th bit
    inline __hostdev__ uint32_t countOn(uint32_t i) const
    {
        uint32_t n = i >> 6, sum = CountOn(mWords[n] & ((uint64_t(1) << (i & 63u)) - 1u));
        for (const uint64_t* w = mWords; n--; ++w)
            sum += CountOn(*w);
        return sum;
    }

    template<bool On>
    class Iterator
    {
    public:
        __hostdev__ Iterator()
            : mPos(Mask::SIZE)
            , mParent(nullptr)
        {
        }
        __hostdev__ Iterator(uint32_t pos, const Mask* parent)
            : mPos(pos)
            , mParent(parent)
        {
        }
        Iterator(const Iterator&) = default;
        Iterator&            operator=(const Iterator&) = default;
        __hostdev__ uint32_t operator*() const { return mPos; }
        __hostdev__ uint32_t pos() const { return mPos; }
        __hostdev__ operator bool() const { return mPos != Mask::SIZE; }
        __hostdev__ Iterator& operator++()
        {
            mPos = mParent->findNext<On>(mPos + 1);
            return *this;
        }
        __hostdev__ Iterator operator++(int)
        {
            auto tmp = *this;
            ++(*this);
            return tmp;
        }

    private:
        uint32_t    mPos;
        const Mask* mParent;
    }; // Member class Iterator

    class DenseIterator
    {
    public:
        __hostdev__ DenseIterator(uint32_t pos = Mask::SIZE)
            : mPos(pos)
        {
        }
        DenseIterator&       operator=(const DenseIterator&) = default;
        __hostdev__ uint32_t operator*() const { return mPos; }
        __hostdev__ uint32_t pos() const { return mPos; }
        __hostdev__ operator bool() const { return mPos != Mask::SIZE; }
        __hostdev__ DenseIterator& operator++()
        {
            ++mPos;
            return *this;
        }
        __hostdev__ DenseIterator operator++(int)
        {
            auto tmp = *this;
            ++mPos;
            return tmp;
        }

    private:
        uint32_t mPos;
    }; // Member class DenseIterator

    using OnIterator = Iterator<true>;
    using OffIterator = Iterator<false>;

    __hostdev__ OnIterator beginOn() const { return OnIterator(this->findFirst<true>(), this); }

    __hostdev__ OffIterator beginOff() const { return OffIterator(this->findFirst<false>(), this); }

    __hostdev__ DenseIterator beginAll() const { return DenseIterator(0); }

    /// @brief Initialize all bits to zero.
    __hostdev__ Mask()
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            mWords[i] = 0;
    }
    __hostdev__ Mask(bool on)
    {
        const uint64_t v = on ? ~uint64_t(0) : uint64_t(0);
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            mWords[i] = v;
    }

    /// @brief Copy constructor
    __hostdev__ Mask(const Mask& other)
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            mWords[i] = other.mWords[i];
    }

    /// @brief Return a pointer to the list of words of the bit mask
    __hostdev__ uint64_t*       words() { return mWords; }
    __hostdev__ const uint64_t* words() const { return mWords; }

    /// @brief Assignment operator that works with openvdb::util::NodeMask
    template<typename MaskT = Mask>
    __hostdev__ typename enable_if<!is_same<MaskT, Mask>::value, Mask&>::type operator=(const MaskT& other)
    {
        static_assert(sizeof(Mask) == sizeof(MaskT), "Mismatching sizeof");
        static_assert(WORD_COUNT == MaskT::WORD_COUNT, "Mismatching word count");
        static_assert(LOG2DIM == MaskT::LOG2DIM, "Mismatching LOG2DIM");
        auto* src = reinterpret_cast<const uint64_t*>(&other);
        for (uint64_t *dst = mWords, *end = dst + WORD_COUNT; dst != end; ++dst)
            *dst = *src++;
        return *this;
    }

    __hostdev__ Mask& operator=(const Mask& other)
    {
        memcpy64(mWords, other.mWords, WORD_COUNT);
        return *this;
    }

    __hostdev__ bool operator==(const Mask& other) const
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i) {
            if (mWords[i] != other.mWords[i])
                return false;
        }
        return true;
    }

    __hostdev__ bool operator!=(const Mask& other) const { return !((*this) == other); }

    /// @brief Return true if the given bit is set.
    __hostdev__ bool isOn(uint32_t n) const { return 0 != (mWords[n >> 6] & (uint64_t(1) << (n & 63))); }

    /// @brief Return true if the given bit is NOT set.
    __hostdev__ bool isOff(uint32_t n) const { return 0 == (mWords[n >> 6] & (uint64_t(1) << (n & 63))); }

    /// @brief Return true if all the bits are set in this Mask.
    __hostdev__ bool isOn() const
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            if (mWords[i] != ~uint64_t(0))
                return false;
        return true;
    }

    /// @brief Return true if none of the bits are set in this Mask.
    __hostdev__ bool isOff() const
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            if (mWords[i] != uint64_t(0))
                return false;
        return true;
    }

    /// @brief Set the specified bit on.
    __hostdev__ void setOn(uint32_t n) { mWords[n >> 6] |= uint64_t(1) << (n & 63); }
    /// @brief Set the specified bit off.
    __hostdev__ void setOff(uint32_t n) { mWords[n >> 6] &= ~(uint64_t(1) << (n & 63)); }

#if defined(__CUDACC__) // the following functions only run on the GPU!
    __device__ inline void setOnAtomic(uint32_t n)
    {
        atomicOr(reinterpret_cast<unsigned long long int*>(this) + (n >> 6), 1ull << (n & 63));
    }
    __device__ inline void setOffAtomic(uint32_t n)
    {
        atomicAnd(reinterpret_cast<unsigned long long int*>(this) + (n >> 6), ~(1ull << (n & 63)));
    }
    __device__ inline void setAtomic(uint32_t n, bool on)
    {
        on ? this->setOnAtomic(n) : this->setOffAtomic(n);
    }
#endif
    /// @brief Set the specified bit on or off.
    __hostdev__ void set(uint32_t n, bool on)
    {
#if 1 // switch between branchless
        auto& word = mWords[n >> 6];
        n &= 63;
        word &= ~(uint64_t(1) << n);
        word |= uint64_t(on) << n;
#else
        on ? this->setOn(n) : this->setOff(n);
#endif
    }

    /// @brief Set all bits on
    __hostdev__ void setOn()
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            mWords[i] = ~uint64_t(0);
    }

    /// @brief Set all bits off
    __hostdev__ void setOff()
    {
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            mWords[i] = uint64_t(0);
    }

    /// @brief Set all bits off
    __hostdev__ void set(bool on)
    {
        const uint64_t v = on ? ~uint64_t(0) : uint64_t(0);
        for (uint32_t i = 0; i < WORD_COUNT; ++i)
            mWords[i] = v;
    }
    /// brief Toggle the state of all bits in the mask
    __hostdev__ void toggle()
    {
        uint32_t n = WORD_COUNT;
        for (auto* w = mWords; n--; ++w)
            *w = ~*w;
    }
    __hostdev__ void toggle(uint32_t n) { mWords[n >> 6] ^= uint64_t(1) << (n & 63); }

    /// @brief Bitwise intersection
    __hostdev__ Mask& operator&=(const Mask& other)
    {
        uint64_t*       w1 = mWords;
        const uint64_t* w2 = other.mWords;
        for (uint32_t n = WORD_COUNT; n--; ++w1, ++w2)
            *w1 &= *w2;
        return *this;
    }
    /// @brief Bitwise union
    __hostdev__ Mask& operator|=(const Mask& other)
    {
        uint64_t*       w1 = mWords;
        const uint64_t* w2 = other.mWords;
        for (uint32_t n = WORD_COUNT; n--; ++w1, ++w2)
            *w1 |= *w2;
        return *this;
    }
    /// @brief Bitwise difference
    __hostdev__ Mask& operator-=(const Mask& other)
    {
        uint64_t*       w1 = mWords;
        const uint64_t* w2 = other.mWords;
        for (uint32_t n = WORD_COUNT; n--; ++w1, ++w2)
            *w1 &= ~*w2;
        return *this;
    }
    /// @brief Bitwise XOR
    __hostdev__ Mask& operator^=(const Mask& other)
    {
        uint64_t*       w1 = mWords;
        const uint64_t* w2 = other.mWords;
        for (uint32_t n = WORD_COUNT; n--; ++w1, ++w2)
            *w1 ^= *w2;
        return *this;
    }

    NANOVDB_HOSTDEV_DISABLE_WARNING
    template<bool ON>
    __hostdev__ uint32_t findFirst() const
    {
        uint32_t        n = 0u;
        const uint64_t* w = mWords;
        for (; n < WORD_COUNT && !(ON ? *w : ~*w); ++w, ++n)
            ;
        return n < WORD_COUNT ? (n << 6) + FindLowestOn(ON ? *w : ~*w) : SIZE;
    }

    NANOVDB_HOSTDEV_DISABLE_WARNING
    template<bool ON>
    __hostdev__ uint32_t findNext(uint32_t start) const
    {
        uint32_t n = start >> 6; // initiate
        if (n >= WORD_COUNT)
            return SIZE; // check for out of bounds
        uint32_t m = start & 63u;
        uint64_t b = ON ? mWords[n] : ~mWords[n];
        if (b & (uint64_t(1u) << m))
            return start; // simple case: start is on/off
        b &= ~uint64_t(0u) << m; // mask out lower bits
        while (!b && ++n < WORD_COUNT)
            b = ON ? mWords[n] : ~mWords[n]; // find next non-zero word
        return b ? (n << 6) + FindLowestOn(b) : SIZE; // catch last word=0
    }

    NANOVDB_HOSTDEV_DISABLE_WARNING
    template<bool ON>
    __hostdev__ uint32_t findPrev(uint32_t start) const
    {
        uint32_t n = start >> 6; // initiate
        if (n >= WORD_COUNT)
            return SIZE; // check for out of bounds
        uint32_t m = start & 63u;
        uint64_t b = ON ? mWords[n] : ~mWords[n];
        if (b & (uint64_t(1u) << m))
            return start; // simple case: start is on/off
        b &= (uint64_t(1u) << m) - 1u; // mask out higher bits
        while (!b && n)
            b = ON ? mWords[--n] : ~mWords[--n]; // find previous non-zero word
        return b ? (n << 6) + FindHighestOn(b) : SIZE; // catch first word=0
    }

private:
    uint64_t mWords[WORD_COUNT];
}; // Mask class

// ----------------------------> Map <--------------------------------------

/// @brief Defines an affine transform and its inverse represented as a 3x3 matrix and a vec3 translation
struct Map
{ // 264B (not 32B aligned!)
    float  mMatF[9]; // 9*4B <- 3x3 matrix
    float  mInvMatF[9]; // 9*4B <- 3x3 matrix
    float  mVecF[3]; // 3*4B <- translation
    float  mTaperF; // 4B, placeholder for taper value
    double mMatD[9]; // 9*8B <- 3x3 matrix
    double mInvMatD[9]; // 9*8B <- 3x3 matrix
    double mVecD[3]; // 3*8B <- translation
    double mTaperD; // 8B, placeholder for taper value

    /// @brief Default constructor for the identity map
    __hostdev__ Map()
        : mMatF{1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f}
        , mInvMatF{1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f}
        , mVecF{0.0f, 0.0f, 0.0f}
        , mTaperF{1.0f}
        , mMatD{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}
        , mInvMatD{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}
        , mVecD{0.0, 0.0, 0.0}
        , mTaperD{1.0}
    {
    }
    __hostdev__ Map(double s, const Vec3d& t = Vec3d(0.0, 0.0, 0.0))
        : mMatF{float(s), 0.0f, 0.0f, 0.0f, float(s), 0.0f, 0.0f, 0.0f, float(s)}
        , mInvMatF{1.0f / float(s), 0.0f, 0.0f, 0.0f, 1.0f / float(s), 0.0f, 0.0f, 0.0f, 1.0f / float(s)}
        , mVecF{float(t[0]), float(t[1]), float(t[2])}
        , mTaperF{1.0f}
        , mMatD{s, 0.0, 0.0, 0.0, s, 0.0, 0.0, 0.0, s}
        , mInvMatD{1.0 / s, 0.0, 0.0, 0.0, 1.0 / s, 0.0, 0.0, 0.0, 1.0 / s}
        , mVecD{t[0], t[1], t[2]}
        , mTaperD{1.0}
    {
    }

    /// @brief Initialize the member data from 3x3 or 4x4 matrices
    /// @note This is not _hostdev__ since then MatT=openvdb::Mat4d will produce warnings
    template<typename MatT, typename Vec3T>
    void set(const MatT& mat, const MatT& invMat, const Vec3T& translate, double taper = 1.0);

    /// @brief Initialize the member data from 4x4 matrices
    /// @note  The last (4th) row of invMat is actually ignored.
    ///        This is not _hostdev__ since then Mat4T=openvdb::Mat4d will produce warnings
    template<typename Mat4T>
    void set(const Mat4T& mat, const Mat4T& invMat, double taper = 1.0) { this->set(mat, invMat, mat[3], taper); }

    template<typename Vec3T>
    void set(double scale, const Vec3T& translation, double taper = 1.0);

    /// @brief Apply the forward affine transformation to a vector using 64bit floating point arithmetics.
    /// @note Typically this operation is used for the scale, rotation and translation of index -> world mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return Forward mapping for affine transformation, i.e. (mat x ijk) + translation
    template<typename Vec3T>
    __hostdev__ Vec3T applyMap(const Vec3T& ijk) const { return matMult(mMatD, mVecD, ijk); }

    /// @brief Apply the forward affine transformation to a vector using 32bit floating point arithmetics.
    /// @note Typically this operation is used for the scale, rotation and translation of index -> world mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return Forward mapping for affine transformation, i.e. (mat x ijk) + translation
    template<typename Vec3T>
    __hostdev__ Vec3T applyMapF(const Vec3T& ijk) const { return matMult(mMatF, mVecF, ijk); }

    /// @brief Apply the linear forward 3x3 transformation to an input 3d vector using 64bit floating point arithmetics,
    ///        e.g. scale and rotation WITHOUT translation.
    /// @note Typically this operation is used for scale and rotation from index -> world mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return linear forward 3x3 mapping of the input vector
    template<typename Vec3T>
    __hostdev__ Vec3T applyJacobian(const Vec3T& ijk) const { return matMult(mMatD, ijk); }

    /// @brief Apply the linear forward 3x3 transformation to an input 3d vector using 32bit floating point arithmetics,
    ///        e.g. scale and rotation WITHOUT translation.
    /// @note Typically this operation is used for scale and rotation from index -> world mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return linear forward 3x3 mapping of the input vector
    template<typename Vec3T>
    __hostdev__ Vec3T applyJacobianF(const Vec3T& ijk) const { return matMult(mMatF, ijk); }

    /// @brief Apply the inverse affine mapping to a vector using 64bit floating point arithmetics.
    /// @note Typically this operation is used for the world -> index mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param xyz 3D vector to be mapped - typically floating point world coordinates
    /// @return Inverse affine mapping of the input @c xyz i.e. (xyz - translation) x mat^-1
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseMap(const Vec3T& xyz) const
    {
        return matMult(mInvMatD, Vec3T(xyz[0] - mVecD[0], xyz[1] - mVecD[1], xyz[2] - mVecD[2]));
    }

    /// @brief Apply the inverse affine mapping to a vector using 32bit floating point arithmetics.
    /// @note Typically this operation is used for the world -> index mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param xyz 3D vector to be mapped - typically floating point world coordinates
    /// @return Inverse affine mapping of the input @c xyz i.e. (xyz - translation) x mat^-1
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseMapF(const Vec3T& xyz) const
    {
        return matMult(mInvMatF, Vec3T(xyz[0] - mVecF[0], xyz[1] - mVecF[1], xyz[2] - mVecF[2]));
    }

    /// @brief Apply the linear inverse 3x3 transformation to an input 3d vector using 64bit floating point arithmetics,
    ///        e.g. inverse scale and inverse rotation WITHOUT translation.
    /// @note Typically this operation is used for scale and rotation from world -> index mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return linear inverse 3x3 mapping of the input vector i.e. xyz x mat^-1
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseJacobian(const Vec3T& xyz) const { return matMult(mInvMatD, xyz); }

    /// @brief Apply the linear inverse 3x3 transformation to an input 3d vector using 32bit floating point arithmetics,
    ///        e.g. inverse scale and inverse rotation WITHOUT translation.
    /// @note Typically this operation is used for scale and rotation from world -> index mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return linear inverse 3x3 mapping of the input vector i.e. xyz x mat^-1
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseJacobianF(const Vec3T& xyz) const { return matMult(mInvMatF, xyz); }

    /// @brief Apply the transposed inverse 3x3 transformation to an input 3d vector using 64bit floating point arithmetics,
    ///        e.g. inverse scale and inverse rotation WITHOUT translation.
    /// @note Typically this operation is used for scale and rotation from world -> index mapping
    /// @tparam Vec3T Template type of the 3D vector to be mapped
    /// @param ijk 3D vector to be mapped - typically floating point index coordinates
    /// @return linear inverse 3x3 mapping of the input vector i.e. xyz x mat^-1
    template<typename Vec3T>
    __hostdev__ Vec3T applyIJT(const Vec3T& xyz) const { return matMultT(mInvMatD, xyz); }
    template<typename Vec3T>
    __hostdev__ Vec3T applyIJTF(const Vec3T& xyz) const { return matMultT(mInvMatF, xyz); }

    /// @brief Return a voxels size in each coordinate direction, measured at the origin
    __hostdev__ Vec3d getVoxelSize() const { return this->applyMap(Vec3d(1)) - this->applyMap(Vec3d(0)); }
}; // Map

template<typename MatT, typename Vec3T>
inline void Map::set(const MatT& mat, const MatT& invMat, const Vec3T& translate, double taper)
{
    float * mf = mMatF, *vf = mVecF, *mif = mInvMatF;
    double *md = mMatD, *vd = mVecD, *mid = mInvMatD;
    mTaperF = static_cast<float>(taper);
    mTaperD = taper;
    for (int i = 0; i < 3; ++i) {
        *vd++ = translate[i]; //translation
        *vf++ = static_cast<float>(translate[i]); //translation
        for (int j = 0; j < 3; ++j) {
            *md++ = mat[j][i]; //transposed
            *mid++ = invMat[j][i];
            *mf++ = static_cast<float>(mat[j][i]); //transposed
            *mif++ = static_cast<float>(invMat[j][i]);
        }
    }
}

template<typename Vec3T>
inline void Map::set(double dx, const Vec3T& trans, double taper)
{
    NANOVDB_ASSERT(dx > 0.0);
    const double mat[3][3] = { {dx, 0.0, 0.0},   // row 0
                               {0.0, dx, 0.0},   // row 1
                               {0.0, 0.0, dx} }; // row 2
    const double idx = 1.0 / dx;
    const double invMat[3][3] = { {idx, 0.0, 0.0},   // row 0
                                  {0.0, idx, 0.0},   // row 1
                                  {0.0, 0.0, idx} }; // row 2
    this->set(mat, invMat, trans, taper);
}

// ----------------------------> GridBlindMetaData <--------------------------------------

struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) GridBlindMetaData
{ // 288 bytes
    static const int      MaxNameSize = 256; // due to NULL termination the maximum length is one less!
    int64_t               mDataOffset; // byte offset to the blind data, relative to this GridBlindMetaData.
    uint64_t              mValueCount; // number of blind values, e.g. point count
    uint32_t              mValueSize;// byte size of each value, e.g. 4 if mDataType=Float and 1 if mDataType=Unknown since that amounts to char
    GridBlindDataSemantic mSemantic; // semantic meaning of the data.
    GridBlindDataClass    mDataClass; // 4 bytes
    GridType              mDataType; // 4 bytes
    char                  mName[MaxNameSize]; // note this includes the NULL termination
    // no padding required for 32 byte alignment

    // disallow copy-construction since methods like blindData and getBlindData uses the this pointer!
    GridBlindMetaData(const GridBlindMetaData&) = delete;

    // disallow copy-assignment since methods like blindData and getBlindData uses the this pointer!
    const GridBlindMetaData& operator=(const GridBlindMetaData&) = delete;

    __hostdev__ void setBlindData(void* blindData) { mDataOffset = PtrDiff(blindData, this); }

    // unsafe
    __hostdev__ const void* blindData() const {return PtrAdd<void>(this, mDataOffset);}

    /// @brief Get a const pointer to the blind data represented by this meta data
    /// @tparam BlindDataT Expected value type of the blind data.
    /// @return Returns NULL if mGridType!=mapToGridType<BlindDataT>(), else a const point of type BlindDataT.
    /// @note Use mDataType=Unknown if BlindDataT is a custom data type unknown to NanoVDB.
    template<typename BlindDataT>
    __hostdev__ const BlindDataT* getBlindData() const
    {
        //if (mDataType != mapToGridType<BlindDataT>()) printf("getBlindData mismatch\n");
        return mDataType == mapToGridType<BlindDataT>() ? PtrAdd<BlindDataT>(this, mDataOffset) : nullptr;
    }

    /// @brief return true if this meta data has a valid combination of semantic, class and value tags
    __hostdev__ bool isValid() const
    {
        auto check = [&]()->bool{
            switch (mDataType){
            case GridType::Unknown: return mValueSize==1u;// i.e. we encode data as mValueCount chars
            case GridType::Float:   return mValueSize==4u;
            case GridType::Double:  return mValueSize==8u;
            case GridType::Int16:   return mValueSize==2u;
            case GridType::Int32:   return mValueSize==4u;
            case GridType::Int64:   return mValueSize==8u;
            case GridType::Vec3f:   return mValueSize==12u;
            case GridType::Vec3d:   return mValueSize==24u;
            case GridType::Half:    return mValueSize==2u;
            case GridType::RGBA8:   return mValueSize==4u;
            case GridType::Fp8:     return mValueSize==1u;
            case GridType::Fp16:    return mValueSize==2u;
            case GridType::Vec4f:   return mValueSize==16u;
            case GridType::Vec4d:   return mValueSize==32u;
            case GridType::Vec3u8:  return mValueSize==3u;
            case GridType::Vec3u16: return mValueSize==6u;
            default: return true;}// all other combinations are valid
        };
        return nanovdb::isValid(mDataClass, mSemantic, mDataType) && check();
    }

    /// @brief return size in bytes of the blind data represented by this blind meta data
    /// @note This size includes possible padding for 32 byte alignment. The actual amount
    ///       of bind data is mValueCount * mValueSize
    __hostdev__ uint64_t blindDataSize() const
    {
        return AlignUp<NANOVDB_DATA_ALIGNMENT>(mValueCount * mValueSize);
    }
}; // GridBlindMetaData

// ----------------------------> NodeTrait <--------------------------------------

/// @brief Struct to derive node type from its level in a given
///        grid, tree or root while preserving constness
template<typename GridOrTreeOrRootT, int LEVEL>
struct NodeTrait;

// Partial template specialization of above Node struct
template<typename GridOrTreeOrRootT>
struct NodeTrait<GridOrTreeOrRootT, 0>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = typename GridOrTreeOrRootT::LeafNodeType;
    using type = typename GridOrTreeOrRootT::LeafNodeType;
};
template<typename GridOrTreeOrRootT>
struct NodeTrait<const GridOrTreeOrRootT, 0>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = const typename GridOrTreeOrRootT::LeafNodeType;
    using type = const typename GridOrTreeOrRootT::LeafNodeType;
};

template<typename GridOrTreeOrRootT>
struct NodeTrait<GridOrTreeOrRootT, 1>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = typename GridOrTreeOrRootT::RootNodeType::ChildNodeType::ChildNodeType;
    using type = typename GridOrTreeOrRootT::RootNodeType::ChildNodeType::ChildNodeType;
};
template<typename GridOrTreeOrRootT>
struct NodeTrait<const GridOrTreeOrRootT, 1>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = const typename GridOrTreeOrRootT::RootNodeType::ChildNodeType::ChildNodeType;
    using type = const typename GridOrTreeOrRootT::RootNodeType::ChildNodeType::ChildNodeType;
};
template<typename GridOrTreeOrRootT>
struct NodeTrait<GridOrTreeOrRootT, 2>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = typename GridOrTreeOrRootT::RootNodeType::ChildNodeType;
    using type = typename GridOrTreeOrRootT::RootNodeType::ChildNodeType;
};
template<typename GridOrTreeOrRootT>
struct NodeTrait<const GridOrTreeOrRootT, 2>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = const typename GridOrTreeOrRootT::RootNodeType::ChildNodeType;
    using type = const typename GridOrTreeOrRootT::RootNodeType::ChildNodeType;
};
template<typename GridOrTreeOrRootT>
struct NodeTrait<GridOrTreeOrRootT, 3>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = typename GridOrTreeOrRootT::RootNodeType;
    using type = typename GridOrTreeOrRootT::RootNodeType;
};

template<typename GridOrTreeOrRootT>
struct NodeTrait<const GridOrTreeOrRootT, 3>
{
    static_assert(GridOrTreeOrRootT::RootNodeType::LEVEL == 3, "Tree depth is not supported");
    using Type = const typename GridOrTreeOrRootT::RootNodeType;
    using type = const typename GridOrTreeOrRootT::RootNodeType;
};

// ----------------------------> Froward decelerations of random access methods <--------------------------------------

template<typename BuildT>
struct GetValue;
template<typename BuildT>
struct SetValue;
template<typename BuildT>
struct SetVoxel;
template<typename BuildT>
struct GetState;
template<typename BuildT>
struct GetDim;
template<typename BuildT>
struct GetLeaf;
template<typename BuildT>
struct ProbeValue;
template<typename BuildT>
struct GetNodeInfo;

// ----------------------------> Grid <--------------------------------------

/*
    The following class and comment is for internal use only

    Memory layout:

    Grid ->       39 x double                          (world bbox and affine transformation)
    Tree -> Root  3 x ValueType + int32_t + N x Tiles  (background,min,max,tileCount + tileCount x Tiles)

    N2 upper InternalNodes each with 2 bit masks, N2 tiles, and min/max values

    N1 lower InternalNodes each with 2 bit masks, N1 tiles, and min/max values

    N0 LeafNodes each with a bit mask, N0 ValueTypes and min/max

    Example layout: ("---" implies it has a custom offset, "..." implies zero or more)
    [GridData][TreeData]---[RootData][ROOT TILES...]---[InternalData<5>]---[InternalData<4>]---[LeafData<3>]---[BLINDMETA...]---[BLIND0]---[BLIND1]---etc.
*/

/// @brief Struct with all the member data of the Grid (useful during serialization of an openvdb grid)
///
/// @note The transform is assumed to be affine (so linear) and have uniform scale! So frustum transforms
///       and non-uniform scaling are not supported (primarily because they complicate ray-tracing in index space)
///
/// @note No client code should (or can) interface with this struct so it can safely be ignored!
struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) GridData
{ // sizeof(GridData) = 672B
    static const int MaxNameSize = 256; // due to NULL termination the maximum length is one less
    uint64_t         mMagic; // 8B (0) magic to validate it is valid grid data.
    uint64_t         mChecksum; // 8B (8). Checksum of grid buffer.
    Version          mVersion; // 4B (16) major, minor, and patch version numbers
    BitFlags<32>     mFlags; // 4B (20). flags for grid.
    uint32_t         mGridIndex; // 4B (24). Index of this grid in the buffer
    uint32_t         mGridCount; // 4B (28). Total number of grids in the buffer
    uint64_t         mGridSize; // 8B (32). byte count of this entire grid occupied in the buffer.
    char             mGridName[MaxNameSize]; // 256B (40)
    Map              mMap; // 264B (296). affine transformation between index and world space in both single and double precision
    BBox<Vec3d>      mWorldBBox; // 48B (560). floating-point AABB of active values in WORLD SPACE (2 x 3 doubles)
    Vec3d            mVoxelSize; // 24B (608). size of a voxel in world units
    GridClass        mGridClass; // 4B (632).
    GridType         mGridType; //  4B (636).
    int64_t          mBlindMetadataOffset; // 8B (640). offset to beginning of GridBlindMetaData structures that follow this grid.
    uint32_t         mBlindMetadataCount; // 4B (648). count of GridBlindMetaData structures that follow this grid.
    uint32_t         mData0; // 4B (652)
    uint64_t         mData1, mData2; // 2x8B (656) padding to 32 B alignment. mData1 is use for the total number of values indexed by an IndexGrid
    /// @brief Use this method to initiate most member dat
    __hostdev__ GridData& operator=(const GridData& other)
    {
        static_assert(8 * 84 == sizeof(GridData), "GridData has unexpected size");
        memcpy64(this, &other, 84);
        return *this;
    }
    __hostdev__ void init(std::initializer_list<GridFlags> list = {GridFlags::IsBreadthFirst},
                          uint64_t                         gridSize = 0u,
                          const Map&                       map = Map(),
                          GridType                         gridType = GridType::Unknown,
                          GridClass                        gridClass = GridClass::Unknown)
    {
#ifdef NANOVDB_USE_NEW_MAGIC_NUMBERS
        mMagic = NANOVDB_MAGIC_GRID;
#else
        mMagic = NANOVDB_MAGIC_NUMBER;
#endif
        mChecksum = ~uint64_t(0);// all 64 bits ON means checksum is disabled
        mVersion = Version();
        mFlags.initMask(list);
        mGridIndex = 0u;
        mGridCount = 1u;
        mGridSize = gridSize;
        mGridName[0] = '\0';
        mMap = map;
        mWorldBBox = BBox<Vec3d>();// invalid bbox
        mVoxelSize = map.getVoxelSize();
        mGridClass = gridClass;
        mGridType = gridType;
        mBlindMetadataOffset = mGridSize; // i.e. no blind data
        mBlindMetadataCount = 0u; // i.e. no blind data
        mData0 = 0u; // zero padding
        mData1 = 0u; // only used for index and point grids
        mData2 = NANOVDB_MAGIC_GRID; // since version 32.6.0 (might be removed in the future)
    }
    /// @brief return true if the magic number and the version are both valid
    __hostdev__ bool isValid() const {
        if (mMagic == NANOVDB_MAGIC_GRID || mData2 == NANOVDB_MAGIC_GRID) return true;
        bool test = mMagic == NANOVDB_MAGIC_NUMBER;// could be GridData or io::FileHeader
        if (test) test = mVersion.isCompatible();
        if (test) test = mGridCount > 0u && mGridIndex < mGridCount;
        if (test) test = mGridClass < GridClass::End && mGridType < GridType::End;
        return test;
    }
    // Set and unset various bit flags
    __hostdev__ void setMinMaxOn(bool on = true) { mFlags.setMask(GridFlags::HasMinMax, on); }
    __hostdev__ void setBBoxOn(bool on = true) { mFlags.setMask(GridFlags::HasBBox, on); }
    __hostdev__ void setLongGridNameOn(bool on = true) { mFlags.setMask(GridFlags::HasLongGridName, on); }
    __hostdev__ void setAverageOn(bool on = true) { mFlags.setMask(GridFlags::HasAverage, on); }
    __hostdev__ void setStdDeviationOn(bool on = true) { mFlags.setMask(GridFlags::HasStdDeviation, on); }
    __hostdev__ bool setGridName(const char* src)
    {
        char *dst = mGridName, *end = dst + MaxNameSize;
        while (*src != '\0' && dst < end - 1)
            *dst++ = *src++;
        while (dst < end)
            *dst++ = '\0';
        return *src == '\0'; // returns true if input grid name is NOT longer than MaxNameSize characters
    }
    // Affine transformations based on double precision
    template<typename Vec3T>
    __hostdev__ Vec3T applyMap(const Vec3T& xyz) const { return mMap.applyMap(xyz); } // Pos: index -> world
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseMap(const Vec3T& xyz) const { return mMap.applyInverseMap(xyz); } // Pos: world -> index
    template<typename Vec3T>
    __hostdev__ Vec3T applyJacobian(const Vec3T& xyz) const { return mMap.applyJacobian(xyz); } // Dir: index -> world
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseJacobian(const Vec3T& xyz) const { return mMap.applyInverseJacobian(xyz); } // Dir: world -> index
    template<typename Vec3T>
    __hostdev__ Vec3T applyIJT(const Vec3T& xyz) const { return mMap.applyIJT(xyz); }
    // Affine transformations based on single precision
    template<typename Vec3T>
    __hostdev__ Vec3T applyMapF(const Vec3T& xyz) const { return mMap.applyMapF(xyz); } // Pos: index -> world
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseMapF(const Vec3T& xyz) const { return mMap.applyInverseMapF(xyz); } // Pos: world -> index
    template<typename Vec3T>
    __hostdev__ Vec3T applyJacobianF(const Vec3T& xyz) const { return mMap.applyJacobianF(xyz); } // Dir: index -> world
    template<typename Vec3T>
    __hostdev__ Vec3T applyInverseJacobianF(const Vec3T& xyz) const { return mMap.applyInverseJacobianF(xyz); } // Dir: world -> index
    template<typename Vec3T>
    __hostdev__ Vec3T applyIJTF(const Vec3T& xyz) const { return mMap.applyIJTF(xyz); }

    // @brief Return a non-const uint8_t pointer to the tree
    __hostdev__ uint8_t* treePtr() { return reinterpret_cast<uint8_t*>(this + 1); }// TreeData is always right after GridData
    //__hostdev__ TreeData* treePtr() { return reinterpret_cast<TreeData*>(this + 1); }// TreeData is always right after GridData

    // @brief Return a const uint8_t pointer to the tree
    __hostdev__ const uint8_t* treePtr() const { return reinterpret_cast<const uint8_t*>(this + 1); }// TreeData is always right after GridData
    //__hostdev__ const TreeData* treePtr() const { return reinterpret_cast<const TreeData*>(this + 1); }// TreeData is always right after GridData

    /// @brief Return a non-const uint8_t pointer to the first node at @c LEVEL
    /// @tparam LEVEL of the node. LEVEL 0 means leaf node and LEVEL 3 means root node
    /// @warning If not nodes exist at @c LEVEL NULL is returned
    template <uint32_t LEVEL>
    __hostdev__ const uint8_t* nodePtr() const
    {
        static_assert(LEVEL >= 0 && LEVEL <= 3, "invalid LEVEL template parameter");
        auto *treeData  = this->treePtr();
        auto nodeOffset = *reinterpret_cast<const uint64_t*>(treeData + 8*LEVEL);// skip LEVEL uint64_t
        return nodeOffset ? PtrAdd<uint8_t>(treeData, nodeOffset) : nullptr;
    }

    /// @brief Return a non-const uint8_t pointer to the first node at @c LEVEL
    /// @tparam LEVEL of the node. LEVEL 0 means leaf node and LEVEL 3 means root node
    /// @warning If not nodes exist at @c LEVEL NULL is returned
    template <uint32_t LEVEL>
    __hostdev__ uint8_t* nodePtr(){return const_cast<uint8_t*>(const_cast<const GridData*>(this)->template nodePtr<LEVEL>());}

    /// @brief Returns a const reference to the blindMetaData at the specified linear offset.
    ///
    /// @warning The linear offset is assumed to be in the valid range
    __hostdev__ const GridBlindMetaData* blindMetaData(uint32_t n) const
    {
        NANOVDB_ASSERT(n < mBlindMetadataCount);
        return PtrAdd<GridBlindMetaData>(this, mBlindMetadataOffset) + n;
    }

    __hostdev__ const char* gridName() const
    {
        if (mFlags.isMaskOn(GridFlags::HasLongGridName)) {// search for first blind meta data that contains a name
            NANOVDB_ASSERT(mBlindMetadataCount > 0);
            for (uint32_t i = 0; i < mBlindMetadataCount; ++i) {
                const auto* metaData = this->blindMetaData(i);// EXTREMELY important to be a pointer
                if (metaData->mDataClass == GridBlindDataClass::GridName) {
                    NANOVDB_ASSERT(metaData->mDataType == GridType::Unknown);
                    return metaData->template getBlindData<const char>();
                }
            }
            NANOVDB_ASSERT(false); // should never hit this!
        }
        return mGridName;
    }

    /// @brief Return memory usage in bytes for this class only.
    __hostdev__ static uint64_t memUsage() { return sizeof(GridData); }

    /// @brief return AABB of active values in world space
    __hostdev__ const BBox<Vec3d>& worldBBox() const { return mWorldBBox; }

    /// @brief return AABB of active values in index space
    __hostdev__ const CoordBBox& indexBBox() const {return *(const CoordBBox*)(this->nodePtr<3>());}

    /// @brief return the root table has size
    __hostdev__ uint32_t rootTableSize() const {
        if (const uint8_t *root = this->nodePtr<3>()) {
            return *(const uint32_t*)(root + sizeof(CoordBBox));
        }
        return 0u;
    }

    /// @brief test if the grid is empty, e.i the root table has size 0
    /// @return  true if this grid contains not data whatsoever
    __hostdev__ bool isEmpty() const {return this->rootTableSize() == 0u;}

    /// @brief  return true if RootData follows TreeData in memory without any extra padding
    /// @details TreeData is always following right after GridData, but the same might not be true for RootData
    __hostdev__ bool isRootConnected() const { return *(const uint64_t*)((const char*)(this + 1) + 24) == 64u;}
}; // GridData

// Forward declaration of accelerated random access class
template<typename BuildT, int LEVEL0 = -1, int LEVEL1 = -1, int LEVEL2 = -1>
class ReadAccessor;

template<typename BuildT>
using DefaultReadAccessor = ReadAccessor<BuildT, 0, 1, 2>;

/// @brief Highest level of the data structure. Contains a tree and a world->index
///        transform (that currently only supports uniform scaling and translation).
///
/// @note This the API of this class to interface with client code
template<typename TreeT>
class Grid : public GridData
{
public:
    using TreeType = TreeT;
    using RootType = typename TreeT::RootType;
    using RootNodeType = RootType;
    using UpperNodeType = typename RootNodeType::ChildNodeType;
    using LowerNodeType = typename UpperNodeType::ChildNodeType;
    using LeafNodeType = typename RootType::LeafNodeType;
    using DataType = GridData;
    using ValueType = typename TreeT::ValueType;
    using BuildType = typename TreeT::BuildType; // in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
    using CoordType = typename TreeT::CoordType;
    using AccessorType = DefaultReadAccessor<BuildType>;

    /// @brief Disallow constructions, copy and assignment
    ///
    /// @note Only a Serializer, defined elsewhere, can instantiate this class
    Grid(const Grid&) = delete;
    Grid& operator=(const Grid&) = delete;
    ~Grid() = delete;

    __hostdev__ Version version() const { return DataType::mVersion; }

    __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }

    __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }

    /// @brief Return memory usage in bytes for this class only.
    //__hostdev__ static uint64_t memUsage() { return sizeof(GridData); }

    /// @brief Return the memory footprint of the entire grid, i.e. including all nodes and blind data
    __hostdev__ uint64_t gridSize() const { return DataType::mGridSize; }

    /// @brief Return index of this grid in the buffer
    __hostdev__ uint32_t gridIndex() const { return DataType::mGridIndex; }

    /// @brief Return total number of grids in the buffer
    __hostdev__ uint32_t gridCount() const { return DataType::mGridCount; }

    /// @brief  @brief Return the total number of values indexed by this IndexGrid
    ///
    /// @note This method is only defined for IndexGrid = NanoGrid<ValueIndex || ValueOnIndex || ValueIndexMask || ValueOnIndexMask>
    template<typename T = BuildType>
    __hostdev__ typename enable_if<BuildTraits<T>::is_index, const uint64_t&>::type
    valueCount() const { return DataType::mData1; }

    /// @brief  @brief Return the total number of points indexed by this PointGrid
    ///
    /// @note This method is only defined for PointGrid = NanoGrid<Point>
    template<typename T = BuildType>
    __hostdev__ typename enable_if<is_same<T, Point>::value, const uint64_t&>::type
    pointCount() const { return DataType::mData1; }

    /// @brief Return a const reference to the tree
    __hostdev__ const TreeT& tree() const { return *reinterpret_cast<const TreeT*>(this->treePtr()); }

    /// @brief Return a non-const reference to the tree
    __hostdev__ TreeT& tree() { return *reinterpret_cast<TreeT*>(this->treePtr()); }

    /// @brief Return a new instance of a ReadAccessor used to access values in this grid
    __hostdev__ AccessorType getAccessor() const { return AccessorType(this->tree().root()); }

    /// @brief Return a const reference to the size of a voxel in world units
    __hostdev__ const Vec3d& voxelSize() const { return DataType::mVoxelSize; }

    /// @brief Return a const reference to the Map for this grid
    __hostdev__ const Map& map() const { return DataType::mMap; }

    /// @brief world to index space transformation
    template<typename Vec3T>
    __hostdev__ Vec3T worldToIndex(const Vec3T& xyz) const { return this->applyInverseMap(xyz); }

    /// @brief index to world space transformation
    template<typename Vec3T>
    __hostdev__ Vec3T indexToWorld(const Vec3T& xyz) const { return this->applyMap(xyz); }

    /// @brief transformation from index space direction to world space direction
    /// @warning assumes dir to be normalized
    template<typename Vec3T>
    __hostdev__ Vec3T indexToWorldDir(const Vec3T& dir) const { return this->applyJacobian(dir); }

    /// @brief transformation from world space direction to index space direction
    /// @warning assumes dir to be normalized
    template<typename Vec3T>
    __hostdev__ Vec3T worldToIndexDir(const Vec3T& dir) const { return this->applyInverseJacobian(dir); }

    /// @brief transform the gradient from index space to world space.
    /// @details Applies the inverse jacobian transform map.
    template<typename Vec3T>
    __hostdev__ Vec3T indexToWorldGrad(const Vec3T& grad) const { return this->applyIJT(grad); }

    /// @brief world to index space transformation
    template<typename Vec3T>
    __hostdev__ Vec3T worldToIndexF(const Vec3T& xyz) const { return this->applyInverseMapF(xyz); }

    /// @brief index to world space transformation
    template<typename Vec3T>
    __hostdev__ Vec3T indexToWorldF(const Vec3T& xyz) const { return this->applyMapF(xyz); }

    /// @brief transformation from index space direction to world space direction
    /// @warning assumes dir to be normalized
    template<typename Vec3T>
    __hostdev__ Vec3T indexToWorldDirF(const Vec3T& dir) const { return this->applyJacobianF(dir); }

    /// @brief transformation from world space direction to index space direction
    /// @warning assumes dir to be normalized
    template<typename Vec3T>
    __hostdev__ Vec3T worldToIndexDirF(const Vec3T& dir) const { return this->applyInverseJacobianF(dir); }

    /// @brief Transforms the gradient from index space to world space.
    /// @details Applies the inverse jacobian transform map.
    template<typename Vec3T>
    __hostdev__ Vec3T indexToWorldGradF(const Vec3T& grad) const { return DataType::applyIJTF(grad); }

    /// @brief Computes a AABB of active values in world space
    //__hostdev__ const BBox<Vec3d>& worldBBox() const { return DataType::mWorldBBox; }

    /// @brief Computes a AABB of active values in index space
    ///
    /// @note This method is returning a floating point bounding box and not a CoordBBox. This makes
    ///       it more useful for clipping rays.
    //__hostdev__ const BBox<CoordType>& indexBBox() const { return this->tree().bbox(); }

    /// @brief Return the total number of active voxels in this tree.
    __hostdev__ uint64_t activeVoxelCount() const { return this->tree().activeVoxelCount(); }

    /// @brief Methods related to the classification of this grid
    __hostdev__ bool             isValid() const { return DataType::isValid(); }
    __hostdev__ const GridType&  gridType() const { return DataType::mGridType; }
    __hostdev__ const GridClass& gridClass() const { return DataType::mGridClass; }
    __hostdev__ bool             isLevelSet() const { return DataType::mGridClass == GridClass::LevelSet; }
    __hostdev__ bool             isFogVolume() const { return DataType::mGridClass == GridClass::FogVolume; }
    __hostdev__ bool             isStaggered() const { return DataType::mGridClass == GridClass::Staggered; }
    __hostdev__ bool             isPointIndex() const { return DataType::mGridClass == GridClass::PointIndex; }
    __hostdev__ bool             isGridIndex() const { return DataType::mGridClass == GridClass::IndexGrid; }
    __hostdev__ bool             isPointData() const { return DataType::mGridClass == GridClass::PointData; }
    __hostdev__ bool             isMask() const { return DataType::mGridClass == GridClass::Topology; }
    __hostdev__ bool             isUnknown() const { return DataType::mGridClass == GridClass::Unknown; }
    __hostdev__ bool             hasMinMax() const { return DataType::mFlags.isMaskOn(GridFlags::HasMinMax); }
    __hostdev__ bool             hasBBox() const { return DataType::mFlags.isMaskOn(GridFlags::HasBBox); }
    __hostdev__ bool             hasLongGridName() const { return DataType::mFlags.isMaskOn(GridFlags::HasLongGridName); }
    __hostdev__ bool             hasAverage() const { return DataType::mFlags.isMaskOn(GridFlags::HasAverage); }
    __hostdev__ bool             hasStdDeviation() const { return DataType::mFlags.isMaskOn(GridFlags::HasStdDeviation); }
    __hostdev__ bool             isBreadthFirst() const { return DataType::mFlags.isMaskOn(GridFlags::IsBreadthFirst); }

    /// @brief return true if the specified node type is layed out breadth-first in memory and has a fixed size.
    ///        This allows for sequential access to the nodes.
    template<typename NodeT>
    __hostdev__ bool isSequential() const { return NodeT::FIXED_SIZE && this->isBreadthFirst(); }

    /// @brief return true if the specified node level is layed out breadth-first in memory and has a fixed size.
    ///        This allows for sequential access to the nodes.
    template<int LEVEL>
    __hostdev__ bool isSequential() const { return NodeTrait<TreeT, LEVEL>::type::FIXED_SIZE && this->isBreadthFirst(); }

   /// @brief return true if nodes at all levels can safely be accessed with simple linear offsets
    __hostdev__ bool isSequential() const { return UpperNodeType::FIXED_SIZE && LowerNodeType::FIXED_SIZE && LeafNodeType::FIXED_SIZE && this->isBreadthFirst(); }

    /// @brief Return a c-string with the name of this grid
    __hostdev__ const char* gridName() const { return DataType::gridName(); }

    /// @brief Return a c-string with the name of this grid, truncated to 255 characters
    __hostdev__ const char* shortGridName() const { return DataType::mGridName; }

    /// @brief Return checksum of the grid buffer.
    __hostdev__ uint64_t checksum() const { return DataType::mChecksum; }

    /// @brief Return true if this grid is empty, i.e. contains no values or nodes.
    //__hostdev__ bool isEmpty() const { return this->tree().isEmpty(); }

    /// @brief Return the count of blind-data encoded in this grid
    __hostdev__ uint32_t blindDataCount() const { return DataType::mBlindMetadataCount; }

    /// @brief Return the index of the first blind data with specified name if found, otherwise -1.
    __hostdev__ int findBlindData(const char* name) const;

    /// @brief Return the index of the first blind data with specified semantic if found, otherwise -1.
    __hostdev__ int findBlindDataForSemantic(GridBlindDataSemantic semantic) const;

    /// @brief Returns a const pointer to the blindData at the specified linear offset.
    ///
    /// @warning Pointer might be NULL and the linear offset is assumed to be in the valid range
    // this method is deprecated !!!!
    __hostdev__ const void* blindData(uint32_t n) const
    {
        printf("\nnanovdb::Grid::blindData is unsafe and hence deprecated! Please use nanovdb::Grid::getBlindData instead.\n\n");
        NANOVDB_ASSERT(n < DataType::mBlindMetadataCount);
        return this->blindMetaData(n).blindData();
    }

    template <typename BlindDataT>
     __hostdev__ const BlindDataT* getBlindData(uint32_t n) const
    {
        if (n >= DataType::mBlindMetadataCount) return nullptr;// index is out of bounds
        return this->blindMetaData(n).template getBlindData<BlindDataT>();// NULL if mismatching BlindDataT
    }

    template <typename BlindDataT>
     __hostdev__ BlindDataT* getBlindData(uint32_t n)
    {
        if (n >= DataType::mBlindMetadataCount) return nullptr;// index is out of bounds
        return const_cast<BlindDataT*>(this->blindMetaData(n).template getBlindData<BlindDataT>());// NULL if mismatching BlindDataT
    }

    __hostdev__ const GridBlindMetaData& blindMetaData(uint32_t n) const { return *DataType::blindMetaData(n); }

private:
    static_assert(sizeof(GridData) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(GridData) is misaligned");
}; // Class Grid

template<typename TreeT>
__hostdev__ int Grid<TreeT>::findBlindDataForSemantic(GridBlindDataSemantic semantic) const
{
    for (uint32_t i = 0, n = this->blindDataCount(); i < n; ++i) {
        if (this->blindMetaData(i).mSemantic == semantic)
            return int(i);
    }
    return -1;
}

template<typename TreeT>
__hostdev__ int Grid<TreeT>::findBlindData(const char* name) const
{
    auto test = [&](int n) {
        const char* str = this->blindMetaData(n).mName;
        for (int i = 0; i < GridBlindMetaData::MaxNameSize; ++i) {
            if (name[i] != str[i])
                return false;
            if (name[i] == '\0' && str[i] == '\0')
                return true;
        }
        return true; // all len characters matched
    };
    for (int i = 0, n = this->blindDataCount(); i < n; ++i)
        if (test(i))
            return i;
    return -1;
}

// ----------------------------> Tree <--------------------------------------

struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) TreeData
{ // sizeof(TreeData) == 64B
    int64_t  mNodeOffset[4];// 32B, byte offset from this tree to first leaf, lower, upper and root node. A zero offset means no node exists
    uint32_t mNodeCount[3]; // 12B, total number of nodes of type: leaf, lower internal, upper internal
    uint32_t mTileCount[3]; // 12B, total number of active tile values at the lower internal, upper internal and root node levels
    uint64_t mVoxelCount; //    8B, total number of active voxels in the root and all its child nodes.
    // No padding since it's always 32B aligned
    __hostdev__ TreeData& operator=(const TreeData& other)
    {
        static_assert(8 * 8 == sizeof(TreeData), "TreeData has unexpected size");
        memcpy64(this, &other, 8);
        return *this;
    }
    __hostdev__ void setRoot(const void* root) {mNodeOffset[3] = root ? PtrDiff(root, this) : 0;}
    __hostdev__ uint8_t* getRoot() { return mNodeOffset[3] ? PtrAdd<uint8_t>(this, mNodeOffset[3]) : nullptr; }
    __hostdev__ const uint8_t* getRoot() const { return mNodeOffset[3] ? PtrAdd<uint8_t>(this, mNodeOffset[3]) : nullptr; }

    template<typename NodeT>
    __hostdev__ void setFirstNode(const NodeT* node) {mNodeOffset[NodeT::LEVEL] = node ? PtrDiff(node, this) : 0;}

    __hostdev__ bool isEmpty() const  {return  mNodeOffset[3] ? *PtrAdd<uint32_t>(this, mNodeOffset[3] + sizeof(BBox<Coord>)) == 0 : true;}

    /// @brief Return the index bounding box of all the active values in this tree, i.e. in all nodes of the tree
    __hostdev__ CoordBBox bbox() const {return  mNodeOffset[3] ? *PtrAdd<CoordBBox>(this, mNodeOffset[3]) : CoordBBox();}

    /// @brief  return true if RootData is layout out immediately after TreeData in memory
    __hostdev__ bool isRootNext() const {return mNodeOffset[3] ? mNodeOffset[3] == sizeof(TreeData) : false; }
};// TreeData

// ----------------------------> GridTree <--------------------------------------

/// @brief defines a tree type from a grid type while preserving constness
template<typename GridT>
struct GridTree
{
    using Type = typename GridT::TreeType;
    using type = typename GridT::TreeType;
};
template<typename GridT>
struct GridTree<const GridT>
{
    using Type = const typename GridT::TreeType;
    using type = const typename GridT::TreeType;
};

// ----------------------------> Tree <--------------------------------------

/// @brief VDB Tree, which is a thin wrapper around a RootNode.
template<typename RootT>
class Tree : public TreeData
{
    static_assert(RootT::LEVEL == 3, "Tree depth is not supported");
    static_assert(RootT::ChildNodeType::LOG2DIM == 5, "Tree configuration is not supported");
    static_assert(RootT::ChildNodeType::ChildNodeType::LOG2DIM == 4, "Tree configuration is not supported");
    static_assert(RootT::LeafNodeType::LOG2DIM == 3, "Tree configuration is not supported");

public:
    using DataType = TreeData;
    using RootType = RootT;
    using RootNodeType = RootT;
    using UpperNodeType = typename RootNodeType::ChildNodeType;
    using LowerNodeType = typename UpperNodeType::ChildNodeType;
    using LeafNodeType = typename RootType::LeafNodeType;
    using ValueType = typename RootT::ValueType;
    using BuildType = typename RootT::BuildType; // in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
    using CoordType = typename RootT::CoordType;
    using AccessorType = DefaultReadAccessor<BuildType>;

    using Node3 = RootT;
    using Node2 = typename RootT::ChildNodeType;
    using Node1 = typename Node2::ChildNodeType;
    using Node0 = LeafNodeType;

    /// @brief This class cannot be constructed or deleted
    Tree() = delete;
    Tree(const Tree&) = delete;
    Tree& operator=(const Tree&) = delete;
    ~Tree() = delete;

    __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }

    __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }

    /// @brief return memory usage in bytes for the class
    __hostdev__ static uint64_t memUsage() { return sizeof(DataType); }

    __hostdev__ RootT& root()
    {
        RootT* ptr = reinterpret_cast<RootT*>(DataType::getRoot());
        NANOVDB_ASSERT(ptr);
        return *ptr;
    }

    __hostdev__ const RootT& root() const
    {
        const RootT* ptr = reinterpret_cast<const RootT*>(DataType::getRoot());
        NANOVDB_ASSERT(ptr);
        return *ptr;
    }

    __hostdev__ AccessorType getAccessor() const { return AccessorType(this->root()); }

    /// @brief Return the value of the given voxel (regardless of state or location in the tree.)
    __hostdev__ ValueType getValue(const CoordType& ijk) const { return this->root().getValue(ijk); }
    __hostdev__ ValueType getValue(int i, int j, int k) const { return this->root().getValue(CoordType(i, j, k)); }

    /// @brief Return the active state of the given voxel (regardless of state or location in the tree.)
    __hostdev__ bool isActive(const CoordType& ijk) const { return this->root().isActive(ijk); }

    /// @brief Return true if this tree is empty, i.e. contains no values or nodes
    //__hostdev__ bool isEmpty() const { return this->root().isEmpty(); }

    /// @brief Combines the previous two methods in a single call
    __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const { return this->root().probeValue(ijk, v); }

    /// @brief Return a const reference to the background value.
    __hostdev__ const ValueType& background() const { return this->root().background(); }

    /// @brief Sets the extrema values of all the active values in this tree, i.e. in all nodes of the tree
    __hostdev__ void extrema(ValueType& min, ValueType& max) const;

    /// @brief Return a const reference to the index bounding box of all the active values in this tree, i.e. in all nodes of the tree
    //__hostdev__ const BBox<CoordType>& bbox() const { return this->root().bbox(); }

    /// @brief Return the total number of active voxels in this tree.
    __hostdev__ uint64_t activeVoxelCount() const { return DataType::mVoxelCount; }

    /// @brief   Return the total number of active tiles at the specified level of the tree.
    ///
    /// @details level = 1,2,3 corresponds to active tile count in lower internal nodes, upper
    ///          internal nodes, and the root level. Note active values at the leaf level are
    ///          referred to as active voxels (see activeVoxelCount defined above).
    __hostdev__ const uint32_t& activeTileCount(uint32_t level) const
    {
        NANOVDB_ASSERT(level > 0 && level <= 3); // 1, 2, or 3
        return DataType::mTileCount[level - 1];
    }

    template<typename NodeT>
    __hostdev__ uint32_t nodeCount() const
    {
        static_assert(NodeT::LEVEL < 3, "Invalid NodeT");
        return DataType::mNodeCount[NodeT::LEVEL];
    }

    __hostdev__ uint32_t nodeCount(int level) const
    {
        NANOVDB_ASSERT(level < 3);
        return DataType::mNodeCount[level];
    }

    __hostdev__ uint32_t totalNodeCount() const
    {
        return DataType::mNodeCount[0] + DataType::mNodeCount[1] + DataType::mNodeCount[2];
    }

    /// @brief return a pointer to the first node of the specified type
    ///
    /// @warning Note it may return NULL if no nodes exist
    template<typename NodeT>
    __hostdev__ NodeT* getFirstNode()
    {
        const int64_t offset = DataType::mNodeOffset[NodeT::LEVEL];
        return offset ? PtrAdd<NodeT>(this, offset) : nullptr;
    }

    /// @brief return a const pointer to the first node of the specified type
    ///
    /// @warning Note it may return NULL if no nodes exist
    template<typename NodeT>
    __hostdev__ const NodeT* getFirstNode() const
    {
        const int64_t offset = DataType::mNodeOffset[NodeT::LEVEL];
        return offset ? PtrAdd<NodeT>(this, offset) : nullptr;
    }

    /// @brief return a pointer to the first node at the specified level
    ///
    /// @warning Note it may return NULL if no nodes exist
    template<int LEVEL>
    __hostdev__ typename NodeTrait<RootT, LEVEL>::type*
    getFirstNode()
    {
        return this->template getFirstNode<typename NodeTrait<RootT, LEVEL>::type>();
    }

    /// @brief return a const pointer to the first node of the specified level
    ///
    /// @warning Note it may return NULL if no nodes exist
    template<int LEVEL>
    __hostdev__ const typename NodeTrait<RootT, LEVEL>::type*
    getFirstNode() const
    {
        return this->template getFirstNode<typename NodeTrait<RootT, LEVEL>::type>();
    }

    /// @brief Template specializations of getFirstNode
    __hostdev__ LeafNodeType*                             getFirstLeaf() { return this->getFirstNode<LeafNodeType>(); }
    __hostdev__ const LeafNodeType*                       getFirstLeaf() const { return this->getFirstNode<LeafNodeType>(); }
    __hostdev__ typename NodeTrait<RootT, 1>::type*       getFirstLower() { return this->getFirstNode<1>(); }
    __hostdev__ const typename NodeTrait<RootT, 1>::type* getFirstLower() const { return this->getFirstNode<1>(); }
    __hostdev__ typename NodeTrait<RootT, 2>::type*       getFirstUpper() { return this->getFirstNode<2>(); }
    __hostdev__ const typename NodeTrait<RootT, 2>::type* getFirstUpper() const { return this->getFirstNode<2>(); }

    template<typename OpT, typename... ArgsT>
    __hostdev__ auto get(const CoordType& ijk, ArgsT&&... args) const
    {
        return this->root().template get<OpT>(ijk, args...);
    }

    template<typename OpT, typename... ArgsT>
    __hostdev__ auto set(const CoordType& ijk, ArgsT&&... args)
    {
        return this->root().template set<OpT>(ijk, args...);
    }

private:
    static_assert(sizeof(DataType) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(TreeData) is misaligned");

}; // Tree class

template<typename RootT>
__hostdev__ void Tree<RootT>::extrema(ValueType& min, ValueType& max) const
{
    min = this->root().minimum();
    max = this->root().maximum();
}

// --------------------------> RootData <------------------------------------

/// @brief Struct with all the member data of the RootNode (useful during serialization of an openvdb RootNode)
///
/// @note No client code should (or can) interface with this struct so it can safely be ignored!
template<typename ChildT>
struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) RootData
{
    using ValueT = typename ChildT::ValueType;
    using BuildT = typename ChildT::BuildType; // in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
    using CoordT = typename ChildT::CoordType;
    using StatsT = typename ChildT::FloatType;
    static constexpr bool FIXED_SIZE = false;

    /// @brief Return a key based on the coordinates of a voxel
#ifdef NANOVDB_USE_SINGLE_ROOT_KEY
    using KeyT = uint64_t;
    template<typename CoordType>
    __hostdev__ static KeyT CoordToKey(const CoordType& ijk)
    {
        static_assert(sizeof(CoordT) == sizeof(CoordType), "Mismatching sizeof");
        static_assert(32 - ChildT::TOTAL <= 21, "Cannot use 64 bit root keys");
        return (KeyT(uint32_t(ijk[2]) >> ChildT::TOTAL)) | //       z is the lower 21 bits
               (KeyT(uint32_t(ijk[1]) >> ChildT::TOTAL) << 21) | // y is the middle 21 bits
               (KeyT(uint32_t(ijk[0]) >> ChildT::TOTAL) << 42); //  x is the upper 21 bits
    }
    __hostdev__ static CoordT KeyToCoord(const KeyT& key)
    {
        static constexpr uint64_t MASK = (1u << 21) - 1; // used to mask out 21 lower bits
        return CoordT(((key >> 42) & MASK) << ChildT::TOTAL, // x are the upper 21 bits
                      ((key >> 21) & MASK) << ChildT::TOTAL, // y are the middle 21 bits
                      (key & MASK) << ChildT::TOTAL); // z are the lower 21 bits
    }
#else
    using KeyT = CoordT;
    __hostdev__ static KeyT   CoordToKey(const CoordT& ijk) { return ijk & ~ChildT::MASK; }
    __hostdev__ static CoordT KeyToCoord(const KeyT& key) { return key; }
#endif
    BBox<CoordT> mBBox; // 24B. AABB of active values in index space.
    uint32_t     mTableSize; // 4B. number of tiles and child pointers in the root node

    ValueT mBackground; // background value, i.e. value of any unset voxel
    ValueT mMinimum; // typically 4B, minimum of all the active values
    ValueT mMaximum; // typically 4B, maximum of all the active values
    StatsT mAverage; // typically 4B, average of all the active values in this node and its child nodes
    StatsT mStdDevi; // typically 4B, standard deviation of all the active values in this node and its child nodes

    /// @brief Return padding of this class in bytes, due to aliasing and 32B alignment
    ///
    /// @note The extra bytes are not necessarily at the end, but can come from aliasing of individual data members.
    __hostdev__ static constexpr uint32_t padding()
    {
        return sizeof(RootData) - (24 + 4 + 3 * sizeof(ValueT) + 2 * sizeof(StatsT));
    }

    struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) Tile
    {
        template<typename CoordType>
        __hostdev__ void setChild(const CoordType& k, const void* ptr, const RootData* data)
        {
            key = CoordToKey(k);
            state = false;
            child = PtrDiff(ptr, data);
        }
        template<typename CoordType, typename ValueType>
        __hostdev__ void setValue(const CoordType& k, bool s, const ValueType& v)
        {
            key = CoordToKey(k);
            state = s;
            value = v;
            child = 0;
        }
        __hostdev__ bool   isChild() const { return child != 0; }
        __hostdev__ bool   isValue() const { return child == 0; }
        __hostdev__ bool   isActive() const { return child == 0 && state; }
        __hostdev__ CoordT origin() const { return KeyToCoord(key); }
        KeyT               key; // NANOVDB_USE_SINGLE_ROOT_KEY ? 8B : 12B
        int64_t            child; // 8B. signed byte offset from this node to the child node.  0 means it is a constant tile, so use value.
        uint32_t           state; // 4B. state of tile value
        ValueT             value; // value of tile (i.e. no child node)
    }; // Tile

    /// @brief Returns a non-const reference to the tile at the specified linear offset.
    ///
    /// @warning The linear offset is assumed to be in the valid range
    __hostdev__ const Tile* tile(uint32_t n) const
    {
        NANOVDB_ASSERT(n < mTableSize);
        return reinterpret_cast<const Tile*>(this + 1) + n;
    }
    __hostdev__ Tile* tile(uint32_t n)
    {
        NANOVDB_ASSERT(n < mTableSize);
        return reinterpret_cast<Tile*>(this + 1) + n;
    }

    __hostdev__ Tile* probeTile(const CoordT& ijk)
    {
#if 1 // switch between linear and binary seach
        const auto key = CoordToKey(ijk);
        for (Tile *p = reinterpret_cast<Tile*>(this + 1), *q = p + mTableSize; p < q; ++p)
            if (p->key == key)
                return p;
        return nullptr;
#else // do not enable binary search if tiles are not guaranteed to be sorted!!!!!!
        int32_t low = 0, high = mTableSize; // low is inclusive and high is exclusive
        while (low != high) {
            int         mid = low + ((high - low) >> 1);
            const Tile* tile = &tiles[mid];
            if (tile->key == key) {
                return tile;
            } else if (tile->key < key) {
                low = mid + 1;
            } else {
                high = mid;
            }
        }
        return nullptr;
#endif
    }

    __hostdev__ inline const Tile* probeTile(const CoordT& ijk) const
    {
        return const_cast<RootData*>(this)->probeTile(ijk);
    }

    /// @brief Returns a const reference to the child node in the specified tile.
    ///
    /// @warning A child node is assumed to exist in the specified tile
    __hostdev__ ChildT* getChild(const Tile* tile)
    {
        NANOVDB_ASSERT(tile->child);
        return PtrAdd<ChildT>(this, tile->child);
    }
    __hostdev__ const ChildT* getChild(const Tile* tile) const
    {
        NANOVDB_ASSERT(tile->child);
        return PtrAdd<ChildT>(this, tile->child);
    }

    __hostdev__ const ValueT& getMin() const { return mMinimum; }
    __hostdev__ const ValueT& getMax() const { return mMaximum; }
    __hostdev__ const StatsT& average() const { return mAverage; }
    __hostdev__ const StatsT& stdDeviation() const { return mStdDevi; }

    __hostdev__ void setMin(const ValueT& v) { mMinimum = v; }
    __hostdev__ void setMax(const ValueT& v) { mMaximum = v; }
    __hostdev__ void setAvg(const StatsT& v) { mAverage = v; }
    __hostdev__ void setDev(const StatsT& v) { mStdDevi = v; }

    /// @brief This class cannot be constructed or deleted
    RootData() = delete;
    RootData(const RootData&) = delete;
    RootData& operator=(const RootData&) = delete;
    ~RootData() = delete;
}; // RootData

// --------------------------> RootNode <------------------------------------

/// @brief Top-most node of the VDB tree structure.
template<typename ChildT>
class RootNode : public RootData<ChildT>
{
public:
    using DataType = RootData<ChildT>;
    using ChildNodeType = ChildT;
    using RootType = RootNode<ChildT>; // this allows RootNode to behave like a Tree
    using RootNodeType = RootType;
    using UpperNodeType = ChildT;
    using LowerNodeType = typename UpperNodeType::ChildNodeType;
    using LeafNodeType = typename ChildT::LeafNodeType;
    using ValueType = typename DataType::ValueT;
    using FloatType = typename DataType::StatsT;
    using BuildType = typename DataType::BuildT; // in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool

    using CoordType = typename ChildT::CoordType;
    using BBoxType = BBox<CoordType>;
    using AccessorType = DefaultReadAccessor<BuildType>;
    using Tile = typename DataType::Tile;
    static constexpr bool FIXED_SIZE = DataType::FIXED_SIZE;

    static constexpr uint32_t LEVEL = 1 + ChildT::LEVEL; // level 0 = leaf

    template<typename RootT>
    class BaseIter
    {
    protected:
        using DataT = typename match_const<DataType, RootT>::type;
        using TileT = typename match_const<Tile, RootT>::type;
        DataT*      mData;
        uint32_t    mPos, mSize;
        __hostdev__ BaseIter(DataT* data = nullptr, uint32_t n = 0)
            : mData(data)
            , mPos(0)
            , mSize(n)
        {
        }

    public:
        __hostdev__ operator bool() const { return mPos < mSize; }
        __hostdev__ uint32_t  pos() const { return mPos; }
        __hostdev__ void      next() { ++mPos; }
        __hostdev__ TileT*    tile() const { return mData->tile(mPos); }
        __hostdev__ CoordType getOrigin() const
        {
            NANOVDB_ASSERT(*this);
            return this->tile()->origin();
        }
        __hostdev__ CoordType getCoord() const
        {
            NANOVDB_ASSERT(*this);
            return this->tile()->origin();
        }
    }; // Member class BaseIter

    template<typename RootT>
    class ChildIter : public BaseIter<RootT>
    {
        static_assert(is_same<typename remove_const<RootT>::type, RootNode>::value, "Invalid RootT");
        using BaseT = BaseIter<RootT>;
        using NodeT = typename match_const<ChildT, RootT>::type;

    public:
        __hostdev__ ChildIter()
            : BaseT()
        {
        }
        __hostdev__ ChildIter(RootT* parent)
            : BaseT(parent->data(), parent->tileCount())
        {
            NANOVDB_ASSERT(BaseT::mData);
            while (*this && !this->tile()->isChild())
                this->next();
        }
        __hostdev__ NodeT& operator*() const
        {
            NANOVDB_ASSERT(*this);
            return *BaseT::mData->getChild(this->tile());
        }
        __hostdev__ NodeT* operator->() const
        {
            NANOVDB_ASSERT(*this);
            return BaseT::mData->getChild(this->tile());
        }
        __hostdev__ ChildIter& operator++()
        {
            NANOVDB_ASSERT(BaseT::mData);
            this->next();
            while (*this && this->tile()->isValue())
                this->next();
            return *this;
        }
        __hostdev__ ChildIter operator++(int)
        {
            auto tmp = *this;
            ++(*this);
            return tmp;
        }
    }; // Member class ChildIter

    using ChildIterator      = ChildIter<RootNode>;
    using ConstChildIterator = ChildIter<const RootNode>;

    __hostdev__ ChildIterator       beginChild() { return ChildIterator(this); }
    __hostdev__ ConstChildIterator cbeginChild() const { return ConstChildIterator(this); }

    template<typename RootT>
    class ValueIter : public BaseIter<RootT>
    {
        using BaseT = BaseIter<RootT>;

    public:
        __hostdev__ ValueIter()
            : BaseT()
        {
        }
        __hostdev__ ValueIter(RootT* parent)
            : BaseT(parent->data(), parent->tileCount())
        {
            NANOVDB_ASSERT(BaseT::mData);
            while (*this && this->tile()->isChild())
                this->next();
        }
        __hostdev__ ValueType operator*() const
        {
            NANOVDB_ASSERT(*this);
            return this->tile()->value;
        }
        __hostdev__ bool isActive() const
        {
            NANOVDB_ASSERT(*this);
            return this->tile()->state;
        }
        __hostdev__ ValueIter& operator++()
        {
            NANOVDB_ASSERT(BaseT::mData);
            this->next();
            while (*this && this->tile()->isChild())
                this->next();
            return *this;
        }
        __hostdev__ ValueIter