OpenVDB  9.0.1
NanoVDB.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
4 /*!
5  \file NanoVDB.h
6 
7  \author Ken Museth
8 
9  \date January 8, 2020
10 
11  \brief Implements a light-weight self-contained VDB data-structure in a
12  single file! In other words, this is a significantly watered-down
13  version of the OpenVDB implementation, with few dependencies - so
14  a one-stop-shop for a minimalistic VDB data structure that run on
15  most platforms!
16 
17  \note It is important to note that NanoVDB (by design) is a read-only
18  sparse GPU (and CPU) friendly data structure intended for applications
19  like rendering and collision detection. As such it obviously lacks
20  a lot of the functionality and features of OpenVDB grids. NanoVDB
21  is essentially a compact linearized (or serialized) representation of
22  an OpenVDB tree with getValue methods only. For best performance use
23  the ReadAccessor::getValue method as opposed to the Tree::getValue
24  method. Note that since a ReadAccessor caches previous access patterns
25  it is by design not thread-safe, so use one instantiation per thread
26  (it is very light-weight). Also, it is not safe to copy accessors between
27  the GPU and CPU! In fact, client code should only interface
28  with the API of the Grid class (all other nodes of the NanoVDB data
29  structure can safely be ignored by most client codes)!
30 
31 
32  \warning NanoVDB grids can only be constructed via tools like openToNanoVDB
33  or the GridBuilder. This explains why none of the grid nodes defined below
34  have public constructors or destructors.
35 
36  \details Please see the following paper for more details on the data structure:
37  K. Museth, “VDB: High-Resolution Sparse Volumes with Dynamic Topology”,
38  ACM Transactions on Graphics 32(3), 2013, which can be found here:
39  http://www.museth.org/Ken/Publications_files/Museth_TOG13.pdf
40 
41 
42  Overview: This file implements the following fundamental class that when combined
43  forms the backbone of the VDB tree data structure:
44 
45  Coord- a signed integer coordinate
46  Vec3 - a 3D vector
47  Vec4 - a 4D vector
48  BBox - a bounding box
49  Mask - a bitmask essential to the non-root tree nodes
50  Map - an affine coordinate transformation
51  Grid - contains a Tree and a map for world<->index transformations. Use
52  this class as the main API with client code!
53  Tree - contains a RootNode and getValue methods that should only be used for debugging
54  RootNode - the top-level node of the VDB data structure
55  InternalNode - the internal nodes of the VDB data structure
56  LeafNode - the lowest level tree nodes that encode voxel values and state
57  ReadAccessor - implements accelerated random access operations
58 
59  Semantics: A VDB data structure encodes values and (binary) states associated with
60  signed integer coordinates. Values encoded at the leaf node level are
61  denoted voxel values, and values associated with other tree nodes are referred
62  to as tile values, which by design cover a larger coordinate index domain.
63 
64 
65  Memory layout:
66 
67  GridData is always at the very beginning of the buffer immediately followed by TreeData!
68  The remaining nodes and blind-data are allowed to be scattered thoughout the buffer,
69  though in practice they are arranged as:
70 
71  GridData: 672 bytes (e.g. magic, checksum, major, flags, index, count, size, name, map, world bbox, voxel size, class, type, offset, count)
72 
73  TreeData: 64 bytes (node counts and byte offsets)
74 
75  ... optional padding ...
76 
77  RootData: size depends on ValueType (index bbox, voxel count, tile count, min/max/avg/standard deviation)
78 
79  Array of: RootData::Tile
80 
81  ... optional padding ...
82 
83  Array of: Upper InternalNodes of size 32^3: bbox, two bit masks, 32768 tile values, and min/max/avg/standard deviation values
84 
85  ... optional padding ...
86 
87  Array of: Lower InternalNodes of size 16^3: bbox, two bit masks, 4096 tile values, and min/max/avg/standard deviation values
88 
89  ... optional padding ...
90 
91  Array of: LeafNodes of size 8^3: bbox, bit masks, 512 voxel values, and min/max/avg/standard deviation values
92 
93 
94  Example layout: ("---" implies it has a custom offset, "..." implies zero or more)
95  [GridData(672B)][TreeData(64B)]---[RootData][N x Root::Tile]---[NodeData<5>]---[ModeData<4>]---[LeafData<3>]---[BLINDMETA...]---[BLIND0]---[BLIND1]---etc.
96 
97 */
98 
99 #ifndef NANOVDB_NANOVDB_H_HAS_BEEN_INCLUDED
100 #define NANOVDB_NANOVDB_H_HAS_BEEN_INCLUDED
101 
102 #define NANOVDB_MAGIC_NUMBER 0x304244566f6e614eUL // "NanoVDB0" in hex - little endian (uint64_t)
103 
104 #define NANOVDB_MAJOR_VERSION_NUMBER 32 // reflects changes to the ABI and hence also the file format
105 #define NANOVDB_MINOR_VERSION_NUMBER 3 // reflects changes to the API but not ABI
106 #define NANOVDB_PATCH_VERSION_NUMBER 3 // reflects changes that does not affect the ABI or API
107 
108 // This replaces a Coord key at the root level with a single uint64_t
109 #define USE_SINGLE_ROOT_KEY
110 
111 // This replaces three levels of Coord keys in the ReadAccessor with one Coord
112 //#define USE_SINGLE_ACCESSOR_KEY
113 
114 #define NANOVDB_FPN_BRANCHLESS
115 
116 #define NANOVDB_DATA_ALIGNMENT 32
117 
118 #if !defined(NANOVDB_ALIGN)
119 #define NANOVDB_ALIGN(n) alignas(n)
120 #endif // !defined(NANOVDB_ALIGN)
121 
122 #ifdef __CUDACC_RTC__
123 
124 typedef signed char int8_t;
125 typedef short int16_t;
126 typedef int int32_t;
127 typedef long long int64_t;
128 typedef unsigned char uint8_t;
129 typedef unsigned int uint32_t;
130 typedef unsigned short uint16_t;
131 typedef unsigned long long uint64_t;
132 
133 #define NANOVDB_ASSERT(x)
134 
135 #define UINT64_C(x) (x ## ULL)
136 
137 #else // __CUDACC_RTC__
138 
139 #include <stdlib.h> // for abs in clang7
140 #include <stdint.h> // for types like int32_t etc
141 #include <stddef.h> // for size_t type
142 #include <cassert> // for assert
143 #include <cstdio> // for sprinf
144 #include <cmath> // for sqrt and fma
145 #include <limits> // for numeric_limits
146 
147 // All asserts can be disabled here, even for debug builds
148 #if 1
149 #define NANOVDB_ASSERT(x) assert(x)
150 #else
151 #define NANOVDB_ASSERT(x)
152 #endif
153 
154 #if defined(NANOVDB_USE_INTRINSICS) && defined(_MSC_VER)
155 #include <intrin.h>
156 #pragma intrinsic(_BitScanReverse)
157 #pragma intrinsic(_BitScanForward)
158 #pragma intrinsic(_BitScanReverse64)
159 #pragma intrinsic(_BitScanForward64)
160 #endif
161 
162 #endif // __CUDACC_RTC__
163 
164 #if defined(__CUDACC__) || defined(__HIP__)
165 // Only define __hostdev__ when using NVIDIA CUDA or HIP compiler
166 #define __hostdev__ __host__ __device__
167 #else
168 #define __hostdev__
169 #endif
170 
171 // The following macro will suppress annoying warnings when nvcc
172 // compiles functions that call (host) intrinsics (which is perfectly valid)
173 #if defined(_MSC_VER) && defined(__CUDACC__)
174 #define NANOVDB_HOSTDEV_DISABLE_WARNING __pragma("hd_warning_disable")
175 #elif defined(__GNUC__) && defined(__CUDACC__)
176 #define NANOVDB_HOSTDEV_DISABLE_WARNING _Pragma("hd_warning_disable")
177 #else
178 #define NANOVDB_HOSTDEV_DISABLE_WARNING
179 #endif
180 
181 // A portable implementation of offsetof - unfortunately it doesn't work with static_assert
182 #define NANOVDB_OFFSETOF(CLASS, MEMBER) ((int)(size_t)((char*)&((CLASS*)0)->MEMBER - (char*)0))
183 
184 namespace nanovdb {
185 
186 // --------------------------> Build types <------------------------------------
187 
188 /// @brief Dummy type for a voxel with a binary mask value, e.g. the active state
189 class ValueMask {};
190 
191 /// @brief Dummy type for a 16 bit floating point values
192 class Half {};
193 
194 /// @brief Dummy type for a 4bit quantization of float point values
195 class Fp4 {};
196 
197 /// @brief Dummy type for a 8bit quantization of float point values
198 class Fp8 {};
199 
200 /// @brief Dummy type for a 16bit quantization of float point values
201 class Fp16 {};
202 
203 /// @brief Dummy type for a variable bit quantization of floating point values
204 class FpN {};
205 
206 // --------------------------> GridType <------------------------------------
207 
208 /// @brief List of types that are currently supported by NanoVDB
209 ///
210 /// @note To expand on this list do:
211 /// 1) Add the new type between Unknown and End in the enum below
212 /// 2) Add the new type to OpenToNanoVDB::processGrid that maps OpenVDB types to GridType
213 /// 3) Verify that the ConvertTrait in NanoToOpenVDB.h works correctly with the new type
214 /// 4) Add the new type to mapToGridType (defined below) that maps NanoVDB types to GridType
215 /// 5) Add the new type to toStr (defined below)
216 enum class GridType : uint32_t { Unknown = 0,
217  Float = 1, // single precision floating point value
218  Double = 2,// double precision floating point value
219  Int16 = 3,// half precision signed integer value
220  Int32 = 4,// single precision signed integer value
221  Int64 = 5,// double precision signed integer value
222  Vec3f = 6,// single precision floating 3D vector
223  Vec3d = 7,// double precision floating 3D vector
224  Mask = 8,// no value, just the active state
225  Half = 9,// half precision floating point value
226  UInt32 = 10,// single precision unsigned integer value
227  Boolean = 11,// boolean value, encoded in bit array
228  RGBA8 = 12,// RGBA packed into 32bit word in reverse-order. R in low bits.
229  Fp4 = 13,// 4bit quantization of float point value
230  Fp8 = 14,// 8bit quantization of float point value
231  Fp16 = 15,// 16bit quantization of float point value
232  FpN = 16,// variable bit quantization of floating point value
233  Vec4f = 17,// single precision floating 4D vector
234  Vec4d = 18,// double precision floating 4D vector
235  End = 19 };
236 
237 #ifndef __CUDACC_RTC__
238 /// @brief Retuns a c-string used to describe a GridType
239 inline const char* toStr(GridType gridType)
240 {
241  static const char * LUT[] = { "?", "float", "double" , "int16", "int32",
242  "int64", "Vec3f", "Vec3d", "Mask", "Half",
243  "uint32", "bool", "RGBA8", "Float4", "Float8",
244  "Float16", "FloatN", "Vec4f", "Vec4d", "End" };
245  static_assert( sizeof(LUT)/sizeof(char*) - 1 == int(GridType::End), "Unexpected size of LUT" );
246  return LUT[static_cast<int>(gridType)];
247 }
248 #endif
249 
250 // --------------------------> GridClass <------------------------------------
251 
252 /// @brief Classes (defined in OpenVDB) that are currently supported by NanoVDB
253 enum class GridClass : uint32_t { Unknown = 0,
254  LevelSet = 1, // narrow band level set, e.g. SDF
255  FogVolume = 2, // fog volume, e.g. density
256  Staggered = 3, // staggered MAC grid, e.g. velocity
257  PointIndex = 4, // point index grid
258  PointData = 5, // point data grid
259  Topology = 6, // grid with active states only (no values)
260  VoxelVolume = 7, // volume of geometric cubes, e.g. minecraft
261  End = 8 };
262 
263 #ifndef __CUDACC_RTC__
264 /// @brief Retuns a c-string used to describe a GridClass
265 inline const char* toStr(GridClass gridClass)
266 {
267  static const char * LUT[] = { "?", "SDF", "FOG" , "MAC", "PNTIDX",
268  "PNTDAT", "TOPO", "VOX", "END" };
269  static_assert( sizeof(LUT)/sizeof(char*) - 1 == int(GridClass::End), "Unexpected size of LUT" );
270  return LUT[static_cast<int>(gridClass)];
271 }
272 #endif
273 
274 // --------------------------> GridFlags <------------------------------------
275 
276 /// @brief Grid flags which indicate what extra information is present in the grid buffer.
277 enum class GridFlags : uint32_t {
278  HasLongGridName = 1 << 0,// grid name is longer than 256 characters
279  HasBBox = 1 << 1,// nodes contain bounding-boxes of active values
280  HasMinMax = 1 << 2,// nodes contain min/max of active values
281  HasAverage = 1 << 3,// nodes contain averages of active values
282  HasStdDeviation = 1 << 4,// nodes contain standard deviations of active values
283  IsBreadthFirst = 1 << 5,// nodes are arranged breadth-first in memory
284  End = 1 << 6,
285 };
286 
287 #ifndef __CUDACC_RTC__
288 /// @brief Retuns a c-string used to describe a GridFlags
289 inline const char* toStr(GridFlags gridFlags)
290 {
291  static const char * LUT[] = { "has long grid name",
292  "has bbox",
293  "has min/max",
294  "has average",
295  "has standard deviation",
296  "is breadth-first",
297  "end" };
298  static_assert( 1 << (sizeof(LUT)/sizeof(char*) - 1) == int(GridFlags::End), "Unexpected size of LUT" );
299  return LUT[static_cast<int>(gridFlags)];
300 }
301 #endif
302 
303 // --------------------------> GridBlindData enums <------------------------------------
304 
305 /// @brief Blind-data Classes that are currently supported by NanoVDB
306 enum class GridBlindDataClass : uint32_t { Unknown = 0,
307  IndexArray = 1,
308  AttributeArray = 2,
309  GridName = 3,
310  End = 4 };
311 
312 /// @brief Blind-data Semantics that are currently understood by NanoVDB
313 enum class GridBlindDataSemantic : uint32_t { Unknown = 0,
314  PointPosition = 1,
315  PointColor = 2,
316  PointNormal = 3,
317  PointRadius = 4,
318  PointVelocity = 5,
319  PointId = 6,
320  End = 7 };
321 
322 // --------------------------> is_same <------------------------------------
323 
324 /// @brief C++11 implementation of std::is_same
325 template<typename T1, typename T2>
326 struct is_same
327 {
328  static constexpr bool value = false;
329 };
330 
331 template<typename T>
332 struct is_same<T, T>
333 {
334  static constexpr bool value = true;
335 };
336 
337 // --------------------------> enable_if <------------------------------------
338 
339 /// @brief C++11 implementation of std::enable_if
340 template <bool, typename T = void>
341 struct enable_if
342 {
343 };
344 
345 template <typename T>
346 struct enable_if<true, T>
347 {
348  using type = T;
349 };
350 
351 // --------------------------> is_floating_point <------------------------------------
352 
353 /// @brief C++11 implementation of std::is_floating_point
354 template<typename T>
356 {
358 };
359 
360 // --------------------------> is_specialization <------------------------------------
361 
362 /// @brief Metafunction used to determine if the first template
363 /// parameter is a specialization of the class template
364 /// given in the second template parameter.
365 ///
366 /// @details is_specialization<Vec3<float>, Vec3>::value == true;
367 template<typename AnyType, template<typename...> class TemplateType>
369 {
370  static const bool value = false;
371 };
372 template<typename... Args, template<typename...> class TemplateType>
373 struct is_specialization<TemplateType<Args...>, TemplateType>
374 {
375  static const bool value = true;
376 };
377 
378 // --------------------------> Value Map <------------------------------------
379 
380 /// @brief Maps one type (e.g. the build types above) to other (actual) types
381 template <typename T>
383 {
384  using Type = T;
385  using type = T;
386 };
387 
388 template<>
390 {
391  using Type = bool;
392  using type = bool;
393 };
394 
395 template<>
397 {
398  using Type = float;
399  using type = float;
400 };
401 
402 template<>
404 {
405  using Type = float;
406  using type = float;
407 };
408 
409 template<>
411 {
412  using Type = float;
413  using type = float;
414 };
415 
416 template<>
418 {
419  using Type = float;
420  using type = float;
421 };
422 
423 template<>
425 {
426  using Type = float;
427  using type = float;
428 };
429 
430 // --------------------------> PtrDiff PtrAdd <------------------------------------
431 
432 template <typename T1, typename T2>
433 __hostdev__ inline static int64_t PtrDiff(const T1* p, const T2* q)
434 {
435  NANOVDB_ASSERT(p && q);
436  return reinterpret_cast<const char*>(p) - reinterpret_cast<const char*>(q);
437 }
438 
439 template <typename DstT, typename SrcT>
440 __hostdev__ inline static DstT* PtrAdd(SrcT *p, int64_t offset)
441 {
442  NANOVDB_ASSERT(p);
443  return reinterpret_cast<DstT*>(reinterpret_cast<char*>(p) + offset);
444 }
445 
446 template <typename DstT, typename SrcT>
447 __hostdev__ inline static const DstT* PtrAdd(const SrcT *p, int64_t offset)
448 {
449  NANOVDB_ASSERT(p);
450  return reinterpret_cast<const DstT*>(reinterpret_cast<const char*>(p) + offset);
451 }
452 // --------------------------> Rgba8 <------------------------------------
453 
454 /// @brief 8-bit red, green, blue, alpha packed into 32 bit unsigned int
455 class Rgba8
456 {
457  union {
458  uint8_t c[4];// 4 color channels of red, green, blue and alpha components.
459  uint32_t packed;// 32 bit packed representation
460  } mData;
461 public:
462  static const int SIZE = 4;
463  using ValueType = uint8_t;
464 
465  Rgba8(const Rgba8&) = default;
466  Rgba8(Rgba8&&) = default;
467  Rgba8& operator=(Rgba8&&) = default;
468  Rgba8& operator=(const Rgba8&) = default;
469  __hostdev__ Rgba8() : mData{0,0,0,0} {static_assert(sizeof(uint32_t) == sizeof(Rgba8),"Unexpected sizeof");}
470  __hostdev__ Rgba8(uint8_t r, uint8_t g, uint8_t b, uint8_t a = 255u) : mData{r, g, b, a} {}
471  explicit __hostdev__ Rgba8(uint8_t v) : Rgba8(v,v,v,v) {}
472  __hostdev__ Rgba8(float r, float g, float b, float a = 1.0f)
473  : mData{(uint8_t(0.5f + r * 255.0f)),// round to nearest
474  (uint8_t(0.5f + g * 255.0f)),// round to nearest
475  (uint8_t(0.5f + b * 255.0f)),// round to nearest
476  (uint8_t(0.5f + a * 255.0f))}// round to nearest
477  {
478  }
479  __hostdev__ bool operator<(const Rgba8& rhs) const { return mData.packed < rhs.mData.packed; }
480  __hostdev__ bool operator==(const Rgba8& rhs) const { return mData.packed == rhs.mData.packed; }
481  __hostdev__ float lengthSqr() const
482  {
483  return 0.0000153787005f*(float(mData.c[0])*mData.c[0] +
484  float(mData.c[1])*mData.c[1] +
485  float(mData.c[2])*mData.c[2]);//1/255^2
486  }
487  __hostdev__ float length() const { return sqrtf(this->lengthSqr() ); }
488  __hostdev__ const uint8_t& operator[](int n) const { return mData.c[n]; }
489  __hostdev__ uint8_t& operator[](int n) { return mData.c[n]; }
490  __hostdev__ const uint32_t& packed() const { return mData.packed; }
491  __hostdev__ uint32_t& packed() { return mData.packed; }
492  __hostdev__ const uint8_t& r() const { return mData.c[0]; }
493  __hostdev__ const uint8_t& g() const { return mData.c[1]; }
494  __hostdev__ const uint8_t& b() const { return mData.c[2]; }
495  __hostdev__ const uint8_t& a() const { return mData.c[3]; }
496  __hostdev__ uint8_t& r() { return mData.c[0]; }
497  __hostdev__ uint8_t& g() { return mData.c[1]; }
498  __hostdev__ uint8_t& b() { return mData.c[2]; }
499  __hostdev__ uint8_t& a() { return mData.c[3]; }
500 };// Rgba8
501 
502 using PackedRGBA8 = Rgba8;// for backwards compatibility
503 
504 // --------------------------> isValue(GridType, GridClass) <------------------------------------
505 
506 /// @brief return true if the GridType maps to a floating point value.
507 __hostdev__ inline bool isFloatingPoint(GridType gridType)
508 {
509  return gridType == GridType::Float ||
510  gridType == GridType::Double ||
511  gridType == GridType::Fp4 ||
512  gridType == GridType::Fp8 ||
513  gridType == GridType::Fp16 ||
514  gridType == GridType::FpN;
515 }
516 
517 // --------------------------> isValue(GridType, GridClass) <------------------------------------
518 
519 /// @brief return true if the combination of GridType and GridClass is valid.
520 __hostdev__ inline bool isValid(GridType gridType, GridClass gridClass)
521 {
522  if (gridClass == GridClass::LevelSet || gridClass == GridClass::FogVolume) {
523  return isFloatingPoint(gridType);
524  } else if (gridClass == GridClass::Staggered) {
525  return gridType == GridType::Vec3f || gridType == GridType::Vec3d ||
526  gridType == GridType::Vec4f || gridType == GridType::Vec4d;
527  } else if (gridClass == GridClass::PointIndex || gridClass == GridClass::PointData) {
528  return gridType == GridType::UInt32;
529  } else if (gridClass == GridClass::VoxelVolume) {
530  return gridType == GridType::RGBA8 || gridType == GridType::Float || gridType == GridType::Double || gridType == GridType::Vec3f || gridType == GridType::Vec3d || gridType == GridType::UInt32;
531  }
532  return gridClass < GridClass::End && gridType < GridType::End;// any valid combination
533 }
534 
535 // ----------------------------> Version class <-------------------------------------
536 
537 /// @brief Bit-compacted representation of all three version numbers
538 ///
539 /// @details major is the top 11 bits, minor is the 11 middle bits and patch is the lower 10 bits
540 class Version
541 {
542  uint32_t mData;// 11 + 11 + 10 bit packing of major + minor + patch
543 public:
544  __hostdev__ Version() : mData( uint32_t(NANOVDB_MAJOR_VERSION_NUMBER) << 21 |
545  uint32_t(NANOVDB_MINOR_VERSION_NUMBER) << 10 |
546  uint32_t(NANOVDB_PATCH_VERSION_NUMBER) )
547  {
548  }
549  __hostdev__ Version(uint32_t major, uint32_t minor, uint32_t patch)
550  : mData( major << 21 | minor << 10 | patch )
551  {
552  NANOVDB_ASSERT(major < (1u << 11));// max value of major is 2047
553  NANOVDB_ASSERT(minor < (1u << 11));// max value of minor is 2047
554  NANOVDB_ASSERT(patch < (1u << 10));// max value of patch is 1023
555  }
556  __hostdev__ bool operator==(const Version &rhs) const {return mData == rhs.mData;}
557  __hostdev__ bool operator< (const Version &rhs) const {return mData < rhs.mData;}
558  __hostdev__ bool operator<=(const Version &rhs) const {return mData <= rhs.mData;}
559  __hostdev__ bool operator> (const Version &rhs) const {return mData > rhs.mData;}
560  __hostdev__ bool operator>=(const Version &rhs) const {return mData >= rhs.mData;}
561  __hostdev__ uint32_t id() const { return mData; }
562  __hostdev__ uint32_t getMajor() const { return (mData >> 21) & ((1u << 11) - 1);}
563  __hostdev__ uint32_t getMinor() const { return (mData >> 10) & ((1u << 11) - 1);}
564  __hostdev__ uint32_t getPatch() const { return mData & ((1u << 10) - 1);}
565 
566 #ifndef __CUDACC_RTC__
567  const char* c_str() const
568  {
569  char *buffer = (char*)malloc(4 + 1 + 4 + 1 + 4 + 1);// xxxx.xxxx.xxxx\n
570  sprintf(buffer, "%d.%d.%d", this->getMajor(), this->getMinor(), this->getPatch());
571  return buffer;
572  }
573 #endif
574 };// Version
575 
576 // ----------------------------> Various math functions <-------------------------------------
577 
578 //@{
579 /// Tolerance for floating-point comparison
580 template<typename T>
581 struct Tolerance;
582 template<>
583 struct Tolerance<float>
584 {
585  __hostdev__ static float value() { return 1e-8f; }
586 };
587 template<>
588 struct Tolerance<double>
589 {
590  __hostdev__ static double value() { return 1e-15; }
591 };
592 //@}
593 
594 //@{
595 /// Delta for small floating-point offsets
596 template<typename T>
597 struct Delta;
598 template<>
599 struct Delta<float>
600 {
601  __hostdev__ static float value() { return 1e-5f; }
602 };
603 template<>
604 struct Delta<double>
605 {
606  __hostdev__ static double value() { return 1e-9; }
607 };
608 //@}
609 
610 //@{
611 /// Maximum floating-point values
612 template<typename T>
613 struct Maximum;
614 #if defined(__CUDA_ARCH__) || defined(__HIP__)
615 template<>
616 struct Maximum<int>
617 {
618  __hostdev__ static int value() { return 2147483647; }
619 };
620 template<>
621 struct Maximum<uint32_t>
622 {
623  __hostdev__ static uint32_t value() { return 4294967295; }
624 };
625 template<>
626 struct Maximum<float>
627 {
628  __hostdev__ static float value() { return 1e+38f; }
629 };
630 template<>
631 struct Maximum<double>
632 {
633  __hostdev__ static double value() { return 1e+308; }
634 };
635 #else
636 template<typename T>
637 struct Maximum
638 {
639  static T value() { return std::numeric_limits<T>::max(); }
640 };
641 #endif
642 //@}
643 
644 template<typename Type>
645 __hostdev__ inline bool isApproxZero(const Type& x)
646 {
647  return !(x > Tolerance<Type>::value()) && !(x < -Tolerance<Type>::value());
648 }
649 
650 template<typename Type>
651 __hostdev__ inline Type Min(Type a, Type b)
652 {
653  return (a < b) ? a : b;
654 }
655 __hostdev__ inline int32_t Min(int32_t a, int32_t b)
656 {
657  return int32_t(fminf(float(a), float(b)));
658 }
659 __hostdev__ inline uint32_t Min(uint32_t a, uint32_t b)
660 {
661  return uint32_t(fminf(float(a), float(b)));
662 }
663 __hostdev__ inline float Min(float a, float b)
664 {
665  return fminf(a, b);
666 }
667 __hostdev__ inline double Min(double a, double b)
668 {
669  return fmin(a, b);
670 }
671 template<typename Type>
672 __hostdev__ inline Type Max(Type a, Type b)
673 {
674  return (a > b) ? a : b;
675 }
676 
677 __hostdev__ inline int32_t Max(int32_t a, int32_t b)
678 {
679  return int32_t(fmaxf(float(a), float(b)));
680 }
681 __hostdev__ inline uint32_t Max(uint32_t a, uint32_t b)
682 {
683  return uint32_t(fmaxf(float(a), float(b)));
684 }
685 __hostdev__ inline float Max(float a, float b)
686 {
687  return fmaxf(a, b);
688 }
689 __hostdev__ inline double Max(double a, double b)
690 {
691  return fmax(a, b);
692 }
693 __hostdev__ inline float Clamp(float x, float a, float b)
694 {
695  return Max(Min(x, b), a);
696 }
697 __hostdev__ inline double Clamp(double x, double a, double b)
698 {
699  return Max(Min(x, b), a);
700 }
701 
702 __hostdev__ inline float Fract(float x)
703 {
704  return x - floorf(x);
705 }
706 __hostdev__ inline double Fract(double x)
707 {
708  return x - floor(x);
709 }
710 
711 __hostdev__ inline int32_t Floor(float x)
712 {
713  return int32_t(floorf(x));
714 }
715 __hostdev__ inline int32_t Floor(double x)
716 {
717  return int32_t(floor(x));
718 }
719 
720 __hostdev__ inline int32_t Ceil(float x)
721 {
722  return int32_t(ceilf(x));
723 }
724 __hostdev__ inline int32_t Ceil(double x)
725 {
726  return int32_t(ceil(x));
727 }
728 
729 template<typename T>
730 __hostdev__ inline T Pow2(T x)
731 {
732  return x * x;
733 }
734 
735 template<typename T>
736 __hostdev__ inline T Pow3(T x)
737 {
738  return x * x * x;
739 }
740 
741 template<typename T>
742 __hostdev__ inline T Pow4(T x)
743 {
744  return Pow2(x * x);
745 }
746 template<typename T>
747 __hostdev__ inline T Abs(T x)
748 {
749  return x < 0 ? -x : x;
750 }
751 
752 template<>
753 __hostdev__ inline float Abs(float x)
754 {
755  return fabs(x);
756 }
757 
758 template<>
759 __hostdev__ inline double Abs(double x)
760 {
761  return fabs(x);
762 }
763 
764 template<>
765 __hostdev__ inline int Abs(int x)
766 {
767  return abs(x);
768 }
769 
770 template<typename CoordT, typename RealT, template<typename> class Vec3T>
771 __hostdev__ inline CoordT Round(const Vec3T<RealT>& xyz);
772 
773 template<typename CoordT, template<typename> class Vec3T>
774 __hostdev__ inline CoordT Round(const Vec3T<float>& xyz)
775 {
776  return CoordT(int32_t(rintf(xyz[0])), int32_t(rintf(xyz[1])), int32_t(rintf(xyz[2])));
777  //return CoordT(int32_t(roundf(xyz[0])), int32_t(roundf(xyz[1])), int32_t(roundf(xyz[2])) );
778  //return CoordT(int32_t(floorf(xyz[0] + 0.5f)), int32_t(floorf(xyz[1] + 0.5f)), int32_t(floorf(xyz[2] + 0.5f)));
779 }
780 
781 template<typename CoordT, template<typename> class Vec3T>
782 __hostdev__ inline CoordT Round(const Vec3T<double>& xyz)
783 {
784  return CoordT(int32_t(floor(xyz[0] + 0.5)), int32_t(floor(xyz[1] + 0.5)), int32_t(floor(xyz[2] + 0.5)));
785 }
786 
787 template<typename CoordT, typename RealT, template<typename> class Vec3T>
788 __hostdev__ inline CoordT RoundDown(const Vec3T<RealT>& xyz)
789 {
790  return CoordT(Floor(xyz[0]), Floor(xyz[1]), Floor(xyz[2]));
791 }
792 
793 //@{
794 /// Return the square root of a floating-point value.
795 __hostdev__ inline float Sqrt(float x)
796 {
797  return sqrtf(x);
798 }
799 __hostdev__ inline double Sqrt(double x)
800 {
801  return sqrt(x);
802 }
803 //@}
804 
805 /// Return the sign of the given value as an integer (either -1, 0 or 1).
806 template <typename T>
807 __hostdev__ inline T Sign(const T &x) { return ((T(0) < x)?T(1):T(0)) - ((x < T(0))?T(1):T(0)); }
808 
809 template<typename Vec3T>
810 __hostdev__ inline int MinIndex(const Vec3T& v)
811 {
812 #if 0
813  static const int hashTable[8] = {2, 1, 9, 1, 2, 9, 0, 0}; //9 are dummy values
814  const int hashKey = ((v[0] < v[1]) << 2) + ((v[0] < v[2]) << 1) + (v[1] < v[2]); // ?*4+?*2+?*1
815  return hashTable[hashKey];
816 #else
817  if (v[0] < v[1] && v[0] < v[2])
818  return 0;
819  if (v[1] < v[2])
820  return 1;
821  else
822  return 2;
823 #endif
824 }
825 
826 template<typename Vec3T>
827 __hostdev__ inline int MaxIndex(const Vec3T& v)
828 {
829 #if 0
830  static const int hashTable[8] = {2, 1, 9, 1, 2, 9, 0, 0}; //9 are dummy values
831  const int hashKey = ((v[0] > v[1]) << 2) + ((v[0] > v[2]) << 1) + (v[1] > v[2]); // ?*4+?*2+?*1
832  return hashTable[hashKey];
833 #else
834  if (v[0] > v[1] && v[0] > v[2])
835  return 0;
836  if (v[1] > v[2])
837  return 1;
838  else
839  return 2;
840 #endif
841 }
842 
843 /// @brief round up byteSize to the nearest wordSize, e.g. to align to machine word: AlignUp<sizeof(size_t)(n)
844 ///
845 /// @details both wordSize and byteSize are in byte units
846 template<uint64_t wordSize>
847 __hostdev__ inline uint64_t AlignUp(uint64_t byteCount)
848 {
849  const uint64_t r = byteCount % wordSize;
850  return r ? byteCount - r + wordSize : byteCount;
851 }
852 
853 // ------------------------------> Coord <--------------------------------------
854 
855 // forward decleration so we can define Coord::asVec3s and Coord::asVec3d
856 template<typename> class Vec3;
857 
858 /// @brief Signed (i, j, k) 32-bit integer coordinate class, similar to openvdb::math::Coord
859 class Coord
860 {
861  int32_t mVec[3]; // private member data - three signed index coordinates
862 public:
863  using ValueType = int32_t;
864  using IndexType = uint32_t;
865 
866  /// @brief Initialize all coordinates to zero.
868  : mVec{0, 0, 0}
869  {
870  }
871 
872  /// @brief Initializes all coordinates to the given signed integer.
874  : mVec{n, n, n}
875  {
876  }
877 
878  /// @brief Initializes coordinate to the given signed integers.
880  : mVec{i, j, k}
881  {
882  }
883 
885  : mVec{ptr[0], ptr[1], ptr[2]}
886  {
887  }
888 
889  __hostdev__ int32_t x() const { return mVec[0]; }
890  __hostdev__ int32_t y() const { return mVec[1]; }
891  __hostdev__ int32_t z() const { return mVec[2]; }
892 
893  __hostdev__ int32_t& x() { return mVec[0]; }
894  __hostdev__ int32_t& y() { return mVec[1]; }
895  __hostdev__ int32_t& z() { return mVec[2]; }
896 
897  __hostdev__ static Coord max() { return Coord(int32_t((1u << 31) - 1)); }
898 
899  __hostdev__ static Coord min() { return Coord(-int32_t((1u << 31) - 1) - 1); }
900 
901  __hostdev__ static size_t memUsage() { return sizeof(Coord); }
902 
903  /// @brief Return a const reference to the given Coord component.
904  /// @warning The argument is assumed to be 0, 1, or 2.
905  __hostdev__ const ValueType& operator[](IndexType i) const { return mVec[i]; }
906 
907  /// @brief Return a non-const reference to the given Coord component.
908  /// @warning The argument is assumed to be 0, 1, or 2.
909  __hostdev__ ValueType& operator[](IndexType i) { return mVec[i]; }
910 
911  /// @brief Assignment operator that works with openvdb::Coord
912  template <typename CoordT>
913  __hostdev__ Coord& operator=(const CoordT &other)
914  {
915  static_assert(sizeof(Coord) == sizeof(CoordT), "Mis-matched sizeof");
916  mVec[0] = other[0];
917  mVec[1] = other[1];
918  mVec[2] = other[2];
919  return *this;
920  }
921 
922  /// @brief Return a new instance with coordinates masked by the given unsigned integer.
923  __hostdev__ Coord operator&(IndexType n) const { return Coord(mVec[0] & n, mVec[1] & n, mVec[2] & n); }
924 
925  // @brief Return a new instance with coordinates left-shifted by the given unsigned integer.
926  __hostdev__ Coord operator<<(IndexType n) const { return Coord(mVec[0] << n, mVec[1] << n, mVec[2] << n); }
927 
928  // @brief Return a new instance with coordinates right-shifted by the given unsigned integer.
929  __hostdev__ Coord operator>>(IndexType n) const { return Coord(mVec[0] >> n, mVec[1] >> n, mVec[2] >> n); }
930 
931  /// @brief Return true if this Coord is lexicographically less than the given Coord.
932  __hostdev__ bool operator<(const Coord& rhs) const
933  {
934  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;
935  }
936 
937  // @brief Return true if the Coord components are identical.
938  __hostdev__ bool operator==(const Coord& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; }
939  __hostdev__ bool operator!=(const Coord& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; }
941  {
942  mVec[0] &= n;
943  mVec[1] &= n;
944  mVec[2] &= n;
945  return *this;
946  }
948  {
949  mVec[0] <<= n;
950  mVec[1] <<= n;
951  mVec[2] <<= n;
952  return *this;
953  }
955  {
956  mVec[0] += n;
957  mVec[1] += n;
958  mVec[2] += n;
959  return *this;
960  }
961  __hostdev__ Coord operator+(const Coord& rhs) const { return Coord(mVec[0] + rhs[0], mVec[1] + rhs[1], mVec[2] + rhs[2]); }
962  __hostdev__ Coord operator-(const Coord& rhs) const { return Coord(mVec[0] - rhs[0], mVec[1] - rhs[1], mVec[2] - rhs[2]); }
964  {
965  mVec[0] += rhs[0];
966  mVec[1] += rhs[1];
967  mVec[2] += rhs[2];
968  return *this;
969  }
971  {
972  mVec[0] -= rhs[0];
973  mVec[1] -= rhs[1];
974  mVec[2] -= rhs[2];
975  return *this;
976  }
977 
978  /// @brief Perform a component-wise minimum with the other Coord.
980  {
981  if (other[0] < mVec[0])
982  mVec[0] = other[0];
983  if (other[1] < mVec[1])
984  mVec[1] = other[1];
985  if (other[2] < mVec[2])
986  mVec[2] = other[2];
987  return *this;
988  }
989 
990  /// @brief Perform a component-wise maximum with the other Coord.
992  {
993  if (other[0] > mVec[0])
994  mVec[0] = other[0];
995  if (other[1] > mVec[1])
996  mVec[1] = other[1];
997  if (other[2] > mVec[2])
998  mVec[2] = other[2];
999  return *this;
1000  }
1001 
1003  {
1004  return Coord(mVec[0] + dx, mVec[1] + dy, mVec[2] + dz);
1005  }
1006 
1007  __hostdev__ Coord offsetBy(ValueType n) const { return this->offsetBy(n, n, n); }
1008 
1009  /// Return true if any of the components of @a a are smaller than the
1010  /// corresponding components of @a b.
1011  __hostdev__ static inline bool lessThan(const Coord& a, const Coord& b)
1012  {
1013  return (a[0] < b[0] || a[1] < b[1] || a[2] < b[2]);
1014  }
1015 
1016  /// @brief Return the largest integer coordinates that are not greater
1017  /// than @a xyz (node centered conversion).
1018  template<typename Vec3T>
1019  __hostdev__ static Coord Floor(const Vec3T& xyz) { return Coord(nanovdb::Floor(xyz[0]), nanovdb::Floor(xyz[1]), nanovdb::Floor(xyz[2])); }
1020 
1021  /// @brief Return a hash key derived from the existing coordinates.
1022  /// @details For details on this hash function please see the VDB paper.
1023  template<int Log2N = 3 + 4 + 5>
1024  __hostdev__ uint32_t hash() const { return ((1 << Log2N) - 1) & (mVec[0] * 73856093 ^ mVec[1] * 19349663 ^ mVec[2] * 83492791); }
1025 
1026  /// @brief Return the octant of this Coord
1027  //__hostdev__ size_t octant() const { return (uint32_t(mVec[0])>>31) | ((uint32_t(mVec[1])>>31)<<1) | ((uint32_t(mVec[2])>>31)<<2); }
1028  __hostdev__ uint8_t octant() const { return (uint8_t(bool(mVec[0] & (1u << 31)))) |
1029  (uint8_t(bool(mVec[1] & (1u << 31))) << 1) |
1030  (uint8_t(bool(mVec[2] & (1u << 31))) << 2); }
1031 
1032  /// @brief Return a single precision floating-point vector of this coordinate
1033  __hostdev__ inline Vec3<float> asVec3s() const;
1034 
1035  /// @brief Return a double precision floating-point vector of this coordinate
1036  __hostdev__ inline Vec3<double> asVec3d() const;
1037 }; // Coord class
1038 
1039 // ----------------------------> Vec3 <--------------------------------------
1040 
1041 /// @brief A simple vector class with three double components, similar to openvdb::math::Vec3
1042 template<typename T>
1043 class Vec3
1044 {
1045  T mVec[3];
1046 
1047 public:
1048  static const int SIZE = 3;
1049  using ValueType = T;
1050  Vec3() = default;
1051  __hostdev__ explicit Vec3(T x)
1052  : mVec{x, x, x}
1053  {
1054  }
1055  __hostdev__ Vec3(T x, T y, T z)
1056  : mVec{x, y, z}
1057  {
1058  }
1059  template<typename T2>
1060  __hostdev__ explicit Vec3(const Vec3<T2>& v)
1061  : mVec{T(v[0]), T(v[1]), T(v[2])}
1062  {
1063  }
1064  __hostdev__ explicit Vec3(const Coord& ijk)
1065  : mVec{T(ijk[0]), T(ijk[1]), T(ijk[2])}
1066  {
1067  }
1068  __hostdev__ bool operator==(const Vec3& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2]; }
1069  __hostdev__ bool operator!=(const Vec3& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2]; }
1070  template<typename Vec3T>
1071  __hostdev__ Vec3& operator=(const Vec3T& rhs)
1072  {
1073  mVec[0] = rhs[0];
1074  mVec[1] = rhs[1];
1075  mVec[2] = rhs[2];
1076  return *this;
1077  }
1078  __hostdev__ const T& operator[](int i) const { return mVec[i]; }
1079  __hostdev__ T& operator[](int i) { return mVec[i]; }
1080  template<typename Vec3T>
1081  __hostdev__ T dot(const Vec3T& v) const { return mVec[0] * v[0] + mVec[1] * v[1] + mVec[2] * v[2]; }
1082  template<typename Vec3T>
1083  __hostdev__ Vec3 cross(const Vec3T& v) const
1084  {
1085  return Vec3(mVec[1] * v[2] - mVec[2] * v[1],
1086  mVec[2] * v[0] - mVec[0] * v[2],
1087  mVec[0] * v[1] - mVec[1] * v[0]);
1088  }
1090  {
1091  return mVec[0] * mVec[0] + mVec[1] * mVec[1] + mVec[2] * mVec[2]; // 5 flops
1092  }
1093  __hostdev__ T length() const { return Sqrt(this->lengthSqr()); }
1094  __hostdev__ Vec3 operator-() const { return Vec3(-mVec[0], -mVec[1], -mVec[2]); }
1095  __hostdev__ Vec3 operator*(const Vec3& v) const { return Vec3(mVec[0] * v[0], mVec[1] * v[1], mVec[2] * v[2]); }
1096  __hostdev__ Vec3 operator/(const Vec3& v) const { return Vec3(mVec[0] / v[0], mVec[1] / v[1], mVec[2] / v[2]); }
1097  __hostdev__ Vec3 operator+(const Vec3& v) const { return Vec3(mVec[0] + v[0], mVec[1] + v[1], mVec[2] + v[2]); }
1098  __hostdev__ Vec3 operator-(const Vec3& v) const { return Vec3(mVec[0] - v[0], mVec[1] - v[1], mVec[2] - v[2]); }
1099  __hostdev__ Vec3 operator*(const T& s) const { return Vec3(s * mVec[0], s * mVec[1], s * mVec[2]); }
1100  __hostdev__ Vec3 operator/(const T& s) const { return (T(1) / s) * (*this); }
1102  {
1103  mVec[0] += v[0];
1104  mVec[1] += v[1];
1105  mVec[2] += v[2];
1106  return *this;
1107  }
1109  {
1110  mVec[0] -= v[0];
1111  mVec[1] -= v[1];
1112  mVec[2] -= v[2];
1113  return *this;
1114  }
1116  {
1117  mVec[0] *= s;
1118  mVec[1] *= s;
1119  mVec[2] *= s;
1120  return *this;
1121  }
1122  __hostdev__ Vec3& operator/=(const T& s) { return (*this) *= T(1) / s; }
1123  __hostdev__ Vec3& normalize() { return (*this) /= this->length(); }
1124  /// @brief Perform a component-wise minimum with the other Coord.
1126  {
1127  if (other[0] < mVec[0])
1128  mVec[0] = other[0];
1129  if (other[1] < mVec[1])
1130  mVec[1] = other[1];
1131  if (other[2] < mVec[2])
1132  mVec[2] = other[2];
1133  return *this;
1134  }
1135 
1136  /// @brief Perform a component-wise maximum with the other Coord.
1138  {
1139  if (other[0] > mVec[0])
1140  mVec[0] = other[0];
1141  if (other[1] > mVec[1])
1142  mVec[1] = other[1];
1143  if (other[2] > mVec[2])
1144  mVec[2] = other[2];
1145  return *this;
1146  }
1147  /// @brief Return the smallest vector component
1149  {
1150  return mVec[0] < mVec[1] ? (mVec[0] < mVec[2] ? mVec[0] : mVec[2]) : (mVec[1] < mVec[2] ? mVec[1] : mVec[2]);
1151  }
1152  /// @brief Return the largest vector component
1154  {
1155  return mVec[0] > mVec[1] ? (mVec[0] > mVec[2] ? mVec[0] : mVec[2]) : (mVec[1] > mVec[2] ? mVec[1] : mVec[2]);
1156  }
1157  __hostdev__ Coord floor() const { return Coord(Floor(mVec[0]), Floor(mVec[1]), Floor(mVec[2])); }
1158  __hostdev__ Coord ceil() const { return Coord(Ceil(mVec[0]), Ceil(mVec[1]), Ceil(mVec[2])); }
1159  __hostdev__ Coord round() const { return Coord(Floor(mVec[0] + 0.5), Floor(mVec[1] + 0.5), Floor(mVec[2] + 0.5)); }
1160 }; // Vec3<T>
1161 
1162 template<typename T1, typename T2>
1163 __hostdev__ inline Vec3<T2> operator*(T1 scalar, const Vec3<T2>& vec)
1164 {
1165  return Vec3<T2>(scalar * vec[0], scalar * vec[1], scalar * vec[2]);
1166 }
1167 template<typename T1, typename T2>
1168 __hostdev__ inline Vec3<T2> operator/(T1 scalar, const Vec3<T2>& vec)
1169 {
1170  return Vec3<T2>(scalar / vec[0], scalar / vec[1], scalar / vec[2]);
1171 }
1172 
1177 
1178 /// @brief Return a single precision floating-point vector of this coordinate
1179 __hostdev__ inline Vec3f Coord::asVec3s() const { return Vec3f(float(mVec[0]), float(mVec[1]), float(mVec[2])); }
1180 
1181 /// @brief Return a double precision floating-point vector of this coordinate
1182 __hostdev__ inline Vec3d Coord::asVec3d() const { return Vec3d(double(mVec[0]), double(mVec[1]), double(mVec[2])); }
1183 
1184 // ----------------------------> Vec4 <--------------------------------------
1185 
1186 /// @brief A simple vector class with three double components, similar to openvdb::math::Vec4
1187 template<typename T>
1188 class Vec4
1189 {
1190  T mVec[4];
1191 
1192 public:
1193  static const int SIZE = 4;
1194  using ValueType = T;
1195  Vec4() = default;
1196  __hostdev__ explicit Vec4(T x)
1197  : mVec{x, x, x, x}
1198  {
1199  }
1200  __hostdev__ Vec4(T x, T y, T z, T w)
1201  : mVec{x, y, z, w}
1202  {
1203  }
1204  template<typename T2>
1205  __hostdev__ explicit Vec4(const Vec4<T2>& v)
1206  : mVec{T(v[0]), T(v[1]), T(v[2]), T(v[3])}
1207  {
1208  }
1209  __hostdev__ bool operator==(const Vec4& rhs) const { return mVec[0] == rhs[0] && mVec[1] == rhs[1] && mVec[2] == rhs[2] && mVec[3] == rhs[3]; }
1210  __hostdev__ bool operator!=(const Vec4& rhs) const { return mVec[0] != rhs[0] || mVec[1] != rhs[1] || mVec[2] != rhs[2] || mVec[3] != rhs[3]; }
1211  template<typename Vec4T>
1212  __hostdev__ Vec4& operator=(const Vec4T& rhs)
1213  {
1214  mVec[0] = rhs[0];
1215  mVec[1] = rhs[1];
1216  mVec[2] = rhs[2];
1217  mVec[3] = rhs[3];
1218  return *this;
1219  }
1220  __hostdev__ const T& operator[](int i) const { return mVec[i]; }
1221  __hostdev__ T& operator[](int i) { return mVec[i]; }
1222  template<typename Vec4T>
1223  __hostdev__ T dot(const Vec4T& v) const { return mVec[0] * v[0] + mVec[1] * v[1] + mVec[2] * v[2] + mVec[3] * v[3]; }
1225  {
1226  return mVec[0] * mVec[0] + mVec[1] * mVec[1] + mVec[2] * mVec[2] + mVec[3] * mVec[3]; // 7 flops
1227  }
1228  __hostdev__ T length() const { return Sqrt(this->lengthSqr()); }
1229  __hostdev__ Vec4 operator-() const { return Vec4(-mVec[0], -mVec[1], -mVec[2], -mVec[3]); }
1230  __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]); }
1231  __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]); }
1232  __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]); }
1233  __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]); }
1234  __hostdev__ Vec4 operator*(const T& s) const { return Vec4(s * mVec[0], s * mVec[1], s * mVec[2], s * mVec[3]); }
1235  __hostdev__ Vec4 operator/(const T& s) const { return (T(1) / s) * (*this); }
1237  {
1238  mVec[0] += v[0];
1239  mVec[1] += v[1];
1240  mVec[2] += v[2];
1241  mVec[3] += v[3];
1242  return *this;
1243  }
1245  {
1246  mVec[0] -= v[0];
1247  mVec[1] -= v[1];
1248  mVec[2] -= v[2];
1249  mVec[3] -= v[3];
1250  return *this;
1251  }
1253  {
1254  mVec[0] *= s;
1255  mVec[1] *= s;
1256  mVec[2] *= s;
1257  mVec[3] *= s;
1258  return *this;
1259  }
1260  __hostdev__ Vec4& operator/=(const T& s) { return (*this) *= T(1) / s; }
1261  __hostdev__ Vec4& normalize() { return (*this) /= this->length(); }
1262  /// @brief Perform a component-wise minimum with the other Coord.
1264  {
1265  if (other[0] < mVec[0])
1266  mVec[0] = other[0];
1267  if (other[1] < mVec[1])
1268  mVec[1] = other[1];
1269  if (other[2] < mVec[2])
1270  mVec[2] = other[2];
1271  if (other[3] < mVec[3])
1272  mVec[3] = other[3];
1273  return *this;
1274  }
1275 
1276  /// @brief Perform a component-wise maximum with the other Coord.
1278  {
1279  if (other[0] > mVec[0])
1280  mVec[0] = other[0];
1281  if (other[1] > mVec[1])
1282  mVec[1] = other[1];
1283  if (other[2] > mVec[2])
1284  mVec[2] = other[2];
1285  if (other[3] > mVec[3])
1286  mVec[3] = other[3];
1287  return *this;
1288  }
1289 }; // Vec4<T>
1290 
1291 template<typename T1, typename T2>
1292 __hostdev__ inline Vec4<T2> operator*(T1 scalar, const Vec4<T2>& vec)
1293 {
1294  return Vec4<T2>(scalar * vec[0], scalar * vec[1], scalar * vec[2], scalar * vec[3]);
1295 }
1296 template<typename T1, typename T2>
1297 __hostdev__ inline Vec4<T2> operator/(T1 scalar, const Vec3<T2>& vec)
1298 {
1299  return Vec4<T2>(scalar / vec[0], scalar / vec[1], scalar / vec[2], scalar / vec[3]);
1300 }
1301 
1306 
1307 // ----------------------------> TensorTraits <--------------------------------------
1308 
1311  is_same<T, Rgba8>::value) ? 1 : 0>
1313 
1314 template<typename T>
1315 struct TensorTraits<T, 0>
1316 {
1317  static const int Rank = 0; // i.e. scalar
1318  static const bool IsScalar = true;
1319  static const bool IsVector = false;
1320  static const int Size = 1;
1321  using ElementType = T;
1322  static T scalar(const T& s) { return s; }
1323 };
1324 
1325 template<typename T>
1326 struct TensorTraits<T, 1>
1327 {
1328  static const int Rank = 1; // i.e. vector
1329  static const bool IsScalar = false;
1330  static const bool IsVector = true;
1331  static const int Size = T::SIZE;
1332  using ElementType = typename T::ValueType;
1333  static ElementType scalar(const T& v) { return v.length(); }
1334 };
1335 
1336 // ----------------------------> FloatTraits <--------------------------------------
1337 
1338 template<typename T, int = sizeof(typename TensorTraits<T>::ElementType)>
1340 {
1341  using FloatType = float;
1342 };
1343 
1344 template<typename T>
1345 struct FloatTraits<T, 8>
1346 {
1347  using FloatType = double;
1348 };
1349 
1350 template<>
1351 struct FloatTraits<bool, 1>
1352 {
1353  using FloatType = bool;
1354 };
1355 
1356 template<>
1358 {
1359  using FloatType = bool;
1360 };
1361 
1362 // ----------------------------> mapping ValueType -> GridType <--------------------------------------
1363 
1364 /// @brief Maps from a templated value type to a GridType enum
1365 template<typename BuildT>
1367 {
1368  if (is_same<BuildT, float>::value) { // resolved at compile-time
1369  return GridType::Float;
1370  } else if (is_same<BuildT, double>::value) {
1371  return GridType::Double;
1372  } else if (is_same<BuildT, int16_t>::value) {
1373  return GridType::Int16;
1374  } else if (is_same<BuildT, int32_t>::value) {
1375  return GridType::Int32;
1376  } else if (is_same<BuildT, int64_t>::value) {
1377  return GridType::Int64;
1378  } else if (is_same<BuildT, Vec3f>::value) {
1379  return GridType::Vec3f;
1380  } else if (is_same<BuildT, Vec3d>::value) {
1381  return GridType::Vec3d;
1382  } else if (is_same<BuildT, uint32_t>::value) {
1383  return GridType::UInt32;
1384  } else if (is_same<BuildT, ValueMask>::value) {
1385  return GridType::Mask;
1386  } else if (is_same<BuildT, bool>::value) {
1387  return GridType::Boolean;
1388  } else if (is_same<BuildT, Rgba8>::value) {
1389  return GridType::RGBA8;
1390  } else if (is_same<BuildT, Fp4>::value) {
1391  return GridType::Fp4;
1392  } else if (is_same<BuildT, Fp8>::value) {
1393  return GridType::Fp8;
1394  } else if (is_same<BuildT, Fp16>::value) {
1395  return GridType::Fp16;
1396  } else if (is_same<BuildT, FpN>::value) {
1397  return GridType::FpN;
1398  } else if (is_same<BuildT, Vec4f>::value) {
1399  return GridType::Vec4f;
1400  } else if (is_same<BuildT, Vec4d>::value) {
1401  return GridType::Vec4d;
1402  }
1403  return GridType::Unknown;
1404 }
1405 
1406 // ----------------------------> matMult <--------------------------------------
1407 
1408 template<typename Vec3T>
1409 __hostdev__ inline Vec3T matMult(const float* mat, const Vec3T& xyz)
1410 {
1411  return Vec3T(fmaf(xyz[0], mat[0], fmaf(xyz[1], mat[1], xyz[2] * mat[2])),
1412  fmaf(xyz[0], mat[3], fmaf(xyz[1], mat[4], xyz[2] * mat[5])),
1413  fmaf(xyz[0], mat[6], fmaf(xyz[1], mat[7], xyz[2] * mat[8]))); // 6 fmaf + 3 mult = 9 flops
1414 }
1415 
1416 template<typename Vec3T>
1417 __hostdev__ inline Vec3T matMult(const double* mat, const Vec3T& xyz)
1418 {
1419  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])),
1420  fma(static_cast<double>(xyz[0]), mat[3], fma(static_cast<double>(xyz[1]), mat[4], static_cast<double>(xyz[2]) * mat[5])),
1421  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
1422 }
1423 
1424 template<typename Vec3T>
1425 __hostdev__ inline Vec3T matMult(const float* mat, const float* vec, const Vec3T& xyz)
1426 {
1427  return Vec3T(fmaf(xyz[0], mat[0], fmaf(xyz[1], mat[1], fmaf(xyz[2], mat[2], vec[0]))),
1428  fmaf(xyz[0], mat[3], fmaf(xyz[1], mat[4], fmaf(xyz[2], mat[5], vec[1]))),
1429  fmaf(xyz[0], mat[6], fmaf(xyz[1], mat[7], fmaf(xyz[2], mat[8], vec[2])))); // 9 fmaf = 9 flops
1430 }
1431 
1432 template<typename Vec3T>
1433 __hostdev__ inline Vec3T matMult(const double* mat, const double* vec, const Vec3T& xyz)
1434 {
1435  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]))),
1436  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]))),
1437  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
1438 }
1439 
1440 // matMultT: Multiply with the transpose:
1441 
1442 template<typename Vec3T>
1443 __hostdev__ inline Vec3T matMultT(const float* mat, const Vec3T& xyz)
1444 {
1445  return Vec3T(fmaf(xyz[0], mat[0], fmaf(xyz[1], mat[3], xyz[2] * mat[6])),
1446  fmaf(xyz[0], mat[1], fmaf(xyz[1], mat[4], xyz[2] * mat[7])),
1447  fmaf(xyz[0], mat[2], fmaf(xyz[1], mat[5], xyz[2] * mat[8]))); // 6 fmaf + 3 mult = 9 flops
1448 }
1449 
1450 template<typename Vec3T>
1451 __hostdev__ inline Vec3T matMultT(const double* mat, const Vec3T& xyz)
1452 {
1453  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])),
1454  fma(static_cast<double>(xyz[0]), mat[1], fma(static_cast<double>(xyz[1]), mat[4], static_cast<double>(xyz[2]) * mat[7])),
1455  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
1456 }
1457 
1458 template<typename Vec3T>
1459 __hostdev__ inline Vec3T matMultT(const float* mat, const float* vec, const Vec3T& xyz)
1460 {
1461  return Vec3T(fmaf(xyz[0], mat[0], fmaf(xyz[1], mat[3], fmaf(xyz[2], mat[6], vec[0]))),
1462  fmaf(xyz[0], mat[1], fmaf(xyz[1], mat[4], fmaf(xyz[2], mat[7], vec[1]))),
1463  fmaf(xyz[0], mat[2], fmaf(xyz[1], mat[5], fmaf(xyz[2], mat[8], vec[2])))); // 9 fmaf = 9 flops
1464 }
1465 
1466 template<typename Vec3T>
1467 __hostdev__ inline Vec3T matMultT(const double* mat, const double* vec, const Vec3T& xyz)
1468 {
1469  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]))),
1470  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]))),
1471  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
1472 }
1473 
1474 // ----------------------------> BBox <-------------------------------------
1475 
1476 // Base-class for static polymorphism (cannot be constructed directly)
1477 template<typename Vec3T>
1478 struct BaseBBox
1479 {
1480  Vec3T mCoord[2];
1481  __hostdev__ bool operator==(const BaseBBox& rhs) const { return mCoord[0] == rhs.mCoord[0] && mCoord[1] == rhs.mCoord[1]; };
1482  __hostdev__ bool operator!=(const BaseBBox& rhs) const { return mCoord[0] != rhs.mCoord[0] || mCoord[1] != rhs.mCoord[1]; };
1483  __hostdev__ const Vec3T& operator[](int i) const { return mCoord[i]; }
1484  __hostdev__ Vec3T& operator[](int i) { return mCoord[i]; }
1485  __hostdev__ Vec3T& min() { return mCoord[0]; }
1486  __hostdev__ Vec3T& max() { return mCoord[1]; }
1487  __hostdev__ const Vec3T& min() const { return mCoord[0]; }
1488  __hostdev__ const Vec3T& max() const { return mCoord[1]; }
1489  __hostdev__ Coord& translate(const Vec3T& xyz)
1490  {
1491  mCoord[0] += xyz;
1492  mCoord[1] += xyz;
1493  return *this;
1494  }
1495  // @brief Expand this bounding box to enclose point (i, j, k).
1496  __hostdev__ BaseBBox& expand(const Vec3T& xyz)
1497  {
1498  mCoord[0].minComponent(xyz);
1499  mCoord[1].maxComponent(xyz);
1500  return *this;
1501  }
1502  //__hostdev__ BaseBBox expandBy(typename Vec3T::ValueType padding) const
1503  //{
1504  // return BaseBBox(mCoord[0].offsetBy(-padding),mCoord[1].offsetBy(padding));
1505  //}
1506  __hostdev__ bool isInside(const Vec3T& xyz)
1507  {
1508  if (xyz[0] < mCoord[0][0] || xyz[1] < mCoord[0][1] || xyz[2] < mCoord[0][2])
1509  return false;
1510  if (xyz[0] > mCoord[1][0] || xyz[1] > mCoord[1][1] || xyz[2] > mCoord[1][2])
1511  return false;
1512  return true;
1513  }
1514 
1515 protected:
1517  __hostdev__ BaseBBox(const Vec3T& min, const Vec3T& max)
1518  : mCoord{min, max}
1519  {
1520  }
1521 }; // BaseBBox
1522 
1524 struct BBox;
1525 
1526 /// @brief Partial template specialization for floating point coordinate types.
1527 ///
1528 /// @note Min is inclusive and max is exclusive. If min = max the dimension of
1529 /// the bounding box is zero and therefore it is also empty.
1530 template<typename Vec3T>
1531 struct BBox<Vec3T, true> : public BaseBBox<Vec3T>
1532 {
1533  using Vec3Type = Vec3T;
1534  using ValueType = typename Vec3T::ValueType;
1535  static_assert(is_floating_point<ValueType>::value, "Expected a floating point coordinate type");
1537  using BaseT::mCoord;
1539  : BaseT(Vec3T( Maximum<typename Vec3T::ValueType>::value()),
1540  Vec3T(-Maximum<typename Vec3T::ValueType>::value()))
1541  {
1542  }
1543  __hostdev__ BBox(const Vec3T& min, const Vec3T& max)
1544  : BaseT(min, max)
1545  {
1546  }
1547  __hostdev__ BBox(const Coord& min, const Coord& max)
1548  : BaseT(Vec3T(ValueType(min[0]), ValueType(min[1]), ValueType(min[2])),
1549  Vec3T(ValueType(max[0] + 1), ValueType(max[1] + 1), ValueType(max[2] + 1)))
1550  {
1551  }
1552  __hostdev__ BBox(const BaseBBox<Coord>& bbox) : BBox(bbox[0], bbox[1]) {}
1553  __hostdev__ bool empty() const { return mCoord[0][0] >= mCoord[1][0] ||
1554  mCoord[0][1] >= mCoord[1][1] ||
1555  mCoord[0][2] >= mCoord[1][2]; }
1556  __hostdev__ Vec3T dim() const { return this->empty() ? Vec3T(0) : this->max() - this->min(); }
1557  __hostdev__ bool isInside(const Vec3T& p) const
1558  {
1559  return p[0] > mCoord[0][0] && p[1] > mCoord[0][1] && p[2] > mCoord[0][2] &&
1560  p[0] < mCoord[1][0] && p[1] < mCoord[1][1] && p[2] < mCoord[1][2];
1561  }
1562 };// BBox<Vec3T, true>
1563 
1564 /// @brief Partial template specialization for integer coordinate types
1565 ///
1566 /// @note Both min and max are INCLUDED in the bbox so dim = max - min + 1. So,
1567 /// if min = max the bounding box contains exactly one point and dim = 1!
1568 template<typename CoordT>
1569 struct BBox<CoordT, false> : public BaseBBox<CoordT>
1570 {
1571  static_assert(is_same<int, typename CoordT::ValueType>::value, "Expected \"int\" coordinate type");
1573  using BaseT::mCoord;
1574  /// @brief Iterator over the domain covered by a BBox
1575  /// @details z is the fastest-moving coordinate.
1576  class Iterator
1577  {
1578  const BBox& mBBox;
1579  CoordT mPos;
1580  public:
1582  : mBBox(b)
1583  , mPos(b.min())
1584  {
1585  }
1587  {
1588  if (mPos[2] < mBBox[1][2]) {// this is the most common case
1589  ++mPos[2];
1590  } else if (mPos[1] < mBBox[1][1]) {
1591  mPos[2] = mBBox[0][2];
1592  ++mPos[1];
1593  } else if (mPos[0] <= mBBox[1][0]) {
1594  mPos[2] = mBBox[0][2];
1595  mPos[1] = mBBox[0][1];
1596  ++mPos[0];
1597  }
1598  return *this;
1599  }
1600  __hostdev__ Iterator operator++(int)
1601  {
1602  auto tmp = *this;
1603  ++(*this);
1604  return tmp;
1605  }
1606  /// @brief Return @c true if the iterator still points to a valid coordinate.
1607  __hostdev__ operator bool() const { return mPos[0] <= mBBox[1][0]; }
1608  __hostdev__ const CoordT& operator*() const { return mPos; }
1609  }; // Iterator
1610  __hostdev__ Iterator begin() const { return Iterator{*this}; }
1612  : BaseT(CoordT::max(), CoordT::min())
1613  {
1614  }
1615  __hostdev__ BBox(const CoordT& min, const CoordT& max)
1616  : BaseT(min, max)
1617  {
1618  }
1619  template<typename SplitT>
1620  __hostdev__ BBox(BBox& other, const SplitT&)
1621  : BaseT(other.mCoord[0], other.mCoord[1])
1622  {
1623  NANOVDB_ASSERT(this->is_divisible());
1624  const int n = MaxIndex(this->dim());
1625  mCoord[1][n] = (mCoord[0][n] + mCoord[1][n]) >> 1;
1626  other.mCoord[0][n] = mCoord[1][n] + 1;
1627  }
1628  __hostdev__ bool is_divisible() const { return mCoord[0][0] < mCoord[1][0] &&
1629  mCoord[0][1] < mCoord[1][1] &&
1630  mCoord[0][2] < mCoord[1][2]; }
1631  /// @brief Return true if this bounding box is empty, i.e. uninitialized
1632  __hostdev__ bool empty() const { return mCoord[0][0] > mCoord[1][0] ||
1633  mCoord[0][1] > mCoord[1][1] ||
1634  mCoord[0][2] > mCoord[1][2]; }
1635  __hostdev__ CoordT dim() const { return this->empty() ? Coord(0) : this->max() - this->min() + Coord(1); }
1636  __hostdev__ uint64_t volume() const { auto d = this->dim(); return uint64_t(d[0])*uint64_t(d[1])*uint64_t(d[2]); }
1637  __hostdev__ bool isInside(const CoordT& p) const { return !(CoordT::lessThan(p, this->min()) || CoordT::lessThan(this->max(), p)); }
1638  __hostdev__ bool isInside(const BBox& b) const
1639  {
1640  return !(CoordT::lessThan(b.min(), this->min()) || CoordT::lessThan(this->max(), b.max()));
1641  }
1642 
1643  /// @warning This converts a CoordBBox into a floating-point bounding box which implies that max += 1 !
1644  template<typename RealT>
1646  {
1647  static_assert(is_floating_point<RealT>::value, "CoordBBox::asReal: Expected a floating point coordinate");
1648  return BBox<Vec3<RealT>>(Vec3<RealT>(RealT(mCoord[0][0]), RealT(mCoord[0][1]), RealT(mCoord[0][2])),
1649  Vec3<RealT>(RealT(mCoord[1][0] + 1), RealT(mCoord[1][1] + 1), RealT(mCoord[1][2] + 1)));
1650  }
1651  /// @brief Return a new instance that is expanded by the specified padding.
1652  __hostdev__ BBox expandBy(typename CoordT::ValueType padding) const
1653  {
1654  return BBox(mCoord[0].offsetBy(-padding), mCoord[1].offsetBy(padding));
1655  }
1656 };// BBox<CoordT, false>
1657 
1660 
1661 // -------------------> Find lowest and highest bit in a word <----------------------------
1662 
1663 /// @brief Returns the index of the lowest, i.e. least significant, on bit in the specified 32 bit word
1664 ///
1665 /// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
1667 __hostdev__ static inline uint32_t FindLowestOn(uint32_t v)
1668 {
1669  NANOVDB_ASSERT(v);
1670 #if defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
1671  unsigned long index;
1672  _BitScanForward(&index, v);
1673  return static_cast<uint32_t>(index);
1674 #elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
1675  return static_cast<uint32_t>(__builtin_ctzl(v));
1676 #else
1677  static const unsigned char DeBruijn[32] = {
1678  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};
1679 // disable unary minus on unsigned warning
1680 #if defined(_MSC_VER) && !defined(__NVCC__)
1681 #pragma warning(push)
1682 #pragma warning(disable : 4146)
1683 #endif
1684  return DeBruijn[uint32_t((v & -v) * 0x077CB531U) >> 27];
1685 #if defined(_MSC_VER) && !defined(__NVCC__)
1686 #pragma warning(pop)
1687 #endif
1688 
1689 #endif
1690 }
1691 
1692 /// @brief Returns the index of the highest, i.e. most significant, on bit in the specified 32 bit word
1693 ///
1694 /// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
1696 __hostdev__ static inline uint32_t FindHighestOn(uint32_t v)
1697 {
1698  NANOVDB_ASSERT(v);
1699 #if defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
1700  unsigned long index;
1701  _BitScanReverse(&index, v);
1702  return static_cast<uint32_t>(index);
1703 #elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
1704  return sizeof(unsigned long) * 8 - 1 - __builtin_clzl(v);
1705 
1706 #else
1707  static const unsigned char DeBruijn[32] = {
1708  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};
1709  v |= v >> 1; // first round down to one less than a power of 2
1710  v |= v >> 2;
1711  v |= v >> 4;
1712  v |= v >> 8;
1713  v |= v >> 16;
1714  return DeBruijn[uint32_t(v * 0x07C4ACDDU) >> 27];
1715 #endif
1716 }
1717 
1718 /// @brief Returns the index of the lowest, i.e. least significant, on bit in the specified 64 bit word
1719 ///
1720 /// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
1722 __hostdev__ static inline uint32_t FindLowestOn(uint64_t v)
1723 {
1724  NANOVDB_ASSERT(v);
1725 #if defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
1726  unsigned long index;
1727  _BitScanForward64(&index, v);
1728  return static_cast<uint32_t>(index);
1729 #elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
1730  return static_cast<uint32_t>(__builtin_ctzll(v));
1731 #else
1732  static const unsigned char DeBruijn[64] = {
1733  0, 1, 2, 53, 3, 7, 54, 27, 4, 38, 41, 8, 34, 55, 48, 28,
1734  62, 5, 39, 46, 44, 42, 22, 9, 24, 35, 59, 56, 49, 18, 29, 11,
1735  63, 52, 6, 26, 37, 40, 33, 47, 61, 45, 43, 21, 23, 58, 17, 10,
1736  51, 25, 36, 32, 60, 20, 57, 16, 50, 31, 19, 15, 30, 14, 13, 12,
1737  };
1738 // disable unary minus on unsigned warning
1739 #if defined(_MSC_VER) && !defined(__NVCC__)
1740 #pragma warning(push)
1741 #pragma warning(disable : 4146)
1742 #endif
1743  return DeBruijn[uint64_t((v & -v) * UINT64_C(0x022FDD63CC95386D)) >> 58];
1744 #if defined(_MSC_VER) && !defined(__NVCC__)
1745 #pragma warning(pop)
1746 #endif
1747 
1748 #endif
1749 }
1750 
1751 /// @brief Returns the index of the highest, i.e. most significant, on bit in the specified 64 bit word
1752 ///
1753 /// @warning Assumes that at least one bit is set in the word, i.e. @a v != uint32_t(0)!
1755 __hostdev__ static inline uint32_t FindHighestOn(uint64_t v)
1756 {
1757  NANOVDB_ASSERT(v);
1758 #if defined(_MSC_VER) && defined(NANOVDB_USE_INTRINSICS)
1759  unsigned long index;
1760  _BitScanReverse64(&index, v);
1761  return static_cast<uint32_t>(index);
1762 #elif (defined(__GNUC__) || defined(__clang__)) && defined(NANOVDB_USE_INTRINSICS)
1763  return sizeof(unsigned long) * 8 - 1 - __builtin_clzll(v);
1764 #else
1765  const uint32_t* p = reinterpret_cast<const uint32_t*>(&v);
1766  return p[1] ? 32u + FindHighestOn(p[1]) : FindHighestOn(p[0]);
1767 #endif
1768 }
1769 
1770 // ----------------------------> CountOn <--------------------------------------
1771 
1772 /// @return Number of bits that are on in the specified 64-bit word
1774 __hostdev__ inline uint32_t CountOn(uint64_t v)
1775 {
1776 #if defined(_MSC_VER) && defined(_M_X64)
1777  v = __popcnt64(v);
1778 #elif (defined(__GNUC__) || defined(__clang__))
1779  v = __builtin_popcountll(v);
1780 #else
1781  // Software Implementation
1782  v = v - ((v >> 1) & uint64_t(0x5555555555555555));
1783  v = (v & uint64_t(0x3333333333333333)) + ((v >> 2) & uint64_t(0x3333333333333333));
1784  v = (((v + (v >> 4)) & uint64_t(0xF0F0F0F0F0F0F0F)) * uint64_t(0x101010101010101)) >> 56;
1785 #endif
1786  return static_cast<uint32_t>(v);
1787 }
1788 
1789 // ----------------------------> Mask <--------------------------------------
1790 
1791 /// @brief Bit-mask to encode active states and facilitate sequential iterators
1792 /// and a fast codec for I/O compression.
1793 template<uint32_t LOG2DIM>
1794 class Mask
1795 {
1796  static constexpr uint32_t SIZE = 1U << (3 * LOG2DIM); // Number of bits in mask
1797  static constexpr uint32_t WORD_COUNT = SIZE >> 6; // Number of 64 bit words
1798  uint64_t mWords[WORD_COUNT];
1799 
1800 public:
1801  /// @brief Return the memory footprint in bytes of this Mask
1802  __hostdev__ static size_t memUsage() { return sizeof(Mask); }
1803 
1804  /// @brief Return the number of bits available in this Mask
1805  __hostdev__ static uint32_t bitCount() { return SIZE; }
1806 
1807  /// @brief Return the number of machine words used by this Mask
1808  __hostdev__ static uint32_t wordCount() { return WORD_COUNT; }
1809 
1810  __hostdev__ uint32_t countOn() const
1811  {
1812  uint32_t sum = 0, n = WORD_COUNT;
1813  for (const uint64_t* w = mWords; n--; ++w)
1814  sum += CountOn(*w);
1815  return sum;
1816  }
1817 
1818  class Iterator
1819  {
1820  public:
1822  : mPos(Mask::SIZE)
1823  , mParent(nullptr)
1824  {
1825  }
1826  __hostdev__ Iterator(uint32_t pos, const Mask* parent)
1827  : mPos(pos)
1828  , mParent(parent)
1829  {
1830  }
1831  Iterator& operator=(const Iterator&) = default;
1832  __hostdev__ uint32_t operator*() const { return mPos; }
1833  __hostdev__ operator bool() const { return mPos != Mask::SIZE; }
1835  {
1836  mPos = mParent->findNextOn(mPos + 1);
1837  return *this;
1838  }
1839 
1840  private:
1841  uint32_t mPos;
1842  const Mask* mParent;
1843  }; // Member class MaskIterator
1844 
1845  /// @brief Initialize all bits to zero.
1847  {
1848  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1849  mWords[i] = 0;
1850  }
1852  {
1853  const uint64_t v = on ? ~uint64_t(0) : uint64_t(0);
1854  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1855  mWords[i] = v;
1856  }
1857 
1858  /// @brief Copy constructor
1859  __hostdev__ Mask(const Mask& other)
1860  {
1861  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1862  mWords[i] = other.mWords[i];
1863  }
1864 
1865  /// @brief Return the <i>n</i>th word of the bit mask, for a word of arbitrary size.
1866  template<typename WordT>
1867  __hostdev__ WordT getWord(int n) const
1868  {
1869  NANOVDB_ASSERT(n * 8 * sizeof(WordT) < SIZE);
1870  return reinterpret_cast<const WordT*>(mWords)[n];
1871  }
1872 
1873  /// @brief Assignment operator that works with openvdb::util::NodeMask
1874  template<typename MaskT>
1875  __hostdev__ Mask& operator=(const MaskT& other)
1876  {
1877  static_assert(sizeof(Mask) == sizeof(MaskT), "Mismatching sizeof");
1878  static_assert(WORD_COUNT == MaskT::WORD_COUNT, "Mismatching word count");
1879  static_assert(LOG2DIM == MaskT::LOG2DIM, "Mismatching LOG2DIM");
1880  auto *src = reinterpret_cast<const uint64_t*>(&other);
1881  uint64_t *dst = mWords;
1882  for (uint32_t i = 0; i < WORD_COUNT; ++i) {
1883  *dst++ = *src++;
1884  }
1885  return *this;
1886  }
1887 
1888  __hostdev__ bool operator==(const Mask& other) const
1889  {
1890  for (uint32_t i = 0; i < WORD_COUNT; ++i) {
1891  if (mWords[i] != other.mWords[i]) return false;
1892  }
1893  return true;
1894  }
1895 
1896  __hostdev__ bool operator!=(const Mask& other) const { return !((*this) == other); }
1897 
1898  __hostdev__ Iterator beginOn() const { return Iterator(this->findFirstOn(), this); }
1899 
1900  /// @brief Return true if the given bit is set.
1901  __hostdev__ bool isOn(uint32_t n) const { return 0 != (mWords[n >> 6] & (uint64_t(1) << (n & 63))); }
1902 
1903  __hostdev__ bool isOn() const
1904  {
1905  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1906  if (mWords[i] != ~uint64_t(0))
1907  return false;
1908  return true;
1909  }
1910 
1911  __hostdev__ bool isOff() const
1912  {
1913  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1914  if (mWords[i] != uint64_t(0))
1915  return false;
1916  return true;
1917  }
1918 
1919  /// @brief Set the given bit on.
1920  __hostdev__ void setOn(uint32_t n) { mWords[n >> 6] |= uint64_t(1) << (n & 63); }
1921  __hostdev__ void setOff(uint32_t n) { mWords[n >> 6] &= ~(uint64_t(1) << (n & 63)); }
1922 
1923  __hostdev__ void set(uint32_t n, bool On)
1924  {
1925 #if 1 // switch between branchless
1926  auto &word = mWords[n >> 6];
1927  n &= 63;
1928  word &= ~(uint64_t(1) << n);
1929  word |= uint64_t(On) << n;
1930 #else
1931  On ? this->setOn(n) : this->setOff(n);
1932 #endif
1933  }
1934 
1935  /// @brief Set all bits on
1937  {
1938  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1939  mWords[i] = ~uint64_t(0);
1940  }
1941 
1942  /// @brief Set all bits off
1944  {
1945  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1946  mWords[i] = uint64_t(0);
1947  }
1948 
1949  /// @brief Set all bits off
1950  __hostdev__ void set(bool on)
1951  {
1952  const uint64_t v = on ? ~uint64_t(0) : uint64_t(0);
1953  for (uint32_t i = 0; i < WORD_COUNT; ++i)
1954  mWords[i] = v;
1955  }
1956  /// brief Toggle the state of all bits in the mask
1958  {
1959  uint32_t n = WORD_COUNT;
1960  for (auto* w = mWords; n--; ++w)
1961  *w = ~*w;
1962  }
1963  __hostdev__ void toggle(uint32_t n) { mWords[n >> 6] ^= uint64_t(1) << (n & 63); }
1964 
1965 private:
1966 
1968  __hostdev__ uint32_t findFirstOn() const
1969  {
1970  uint32_t n = 0;
1971  const uint64_t* w = mWords;
1972  for (; n < WORD_COUNT && !*w; ++w, ++n)
1973  ;
1974  return n == WORD_COUNT ? SIZE : (n << 6) + FindLowestOn(*w);
1975  }
1976 
1978  __hostdev__ uint32_t findNextOn(uint32_t start) const
1979  {
1980  uint32_t n = start >> 6; // initiate
1981  if (n >= WORD_COUNT)
1982  return SIZE; // check for out of bounds
1983  uint32_t m = start & 63;
1984  uint64_t b = mWords[n];
1985  if (b & (uint64_t(1) << m))
1986  return start; // simple case: start is on
1987  b &= ~uint64_t(0) << m; // mask out lower bits
1988  while (!b && ++n < WORD_COUNT)
1989  b = mWords[n]; // find next non-zero word
1990  return (!b ? SIZE : (n << 6) + FindLowestOn(b)); // catch last word=0
1991  }
1992 }; // Mask class
1993 
1994 // ----------------------------> Map <--------------------------------------
1995 
1996 /// @brief Defines an affine transform and its inverse represented as a 3x3 matrix and a vec3 translation
1997 struct Map
1998 {
1999  float mMatF[9]; // 9*4B <- 3x3 matrix
2000  float mInvMatF[9]; // 9*4B <- 3x3 matrix
2001  float mVecF[3]; // 3*4B <- translation
2002  float mTaperF; // 4B, placeholder for taper value
2003  double mMatD[9]; // 9*8B <- 3x3 matrix
2004  double mInvMatD[9]; // 9*8B <- 3x3 matrix
2005  double mVecD[3]; // 3*8B <- translation
2006  double mTaperD; // 8B, placeholder for taper value
2007 
2008  // This method can only be called on the host to initialize the member data
2009  template<typename Mat4T>
2010  __hostdev__ void set(const Mat4T& mat, const Mat4T& invMat, double taper);
2011 
2012  template<typename Vec3T>
2013  __hostdev__ Vec3T applyMap(const Vec3T& xyz) const { return matMult(mMatD, mVecD, xyz); }
2014  template<typename Vec3T>
2015  __hostdev__ Vec3T applyMapF(const Vec3T& xyz) const { return matMult(mMatF, mVecF, xyz); }
2016 
2017  template<typename Vec3T>
2018  __hostdev__ Vec3T applyJacobian(const Vec3T& xyz) const { return matMult(mMatD, xyz); }
2019  template<typename Vec3T>
2020  __hostdev__ Vec3T applyJacobianF(const Vec3T& xyz) const { return matMult(mMatF, xyz); }
2021 
2022  template<typename Vec3T>
2023  __hostdev__ Vec3T applyInverseMap(const Vec3T& xyz) const
2024  {
2025  return matMult(mInvMatD, Vec3T(xyz[0] - mVecD[0], xyz[1] - mVecD[1], xyz[2] - mVecD[2]));
2026  }
2027  template<typename Vec3T>
2028  __hostdev__ Vec3T applyInverseMapF(const Vec3T& xyz) const
2029  {
2030  return matMult(mInvMatF, Vec3T(xyz[0] - mVecF[0], xyz[1] - mVecF[1], xyz[2] - mVecF[2]));
2031  }
2032 
2033  template<typename Vec3T>
2034  __hostdev__ Vec3T applyInverseJacobian(const Vec3T& xyz) const { return matMult(mInvMatD, xyz); }
2035  template<typename Vec3T>
2036  __hostdev__ Vec3T applyInverseJacobianF(const Vec3T& xyz) const { return matMult(mInvMatF, xyz); }
2037 
2038  template<typename Vec3T>
2039  __hostdev__ Vec3T applyIJT(const Vec3T& xyz) const { return matMultT(mInvMatD, xyz); }
2040  template<typename Vec3T>
2041  __hostdev__ Vec3T applyIJTF(const Vec3T& xyz) const { return matMultT(mInvMatF, xyz); }
2042 }; // Map
2043 
2044 template<typename Mat4T>
2045 __hostdev__ void Map::set(const Mat4T& mat, const Mat4T& invMat, double taper)
2046 {
2047  float * mf = mMatF, *vf = mVecF;
2048  float* mif = mInvMatF;
2049  double *md = mMatD, *vd = mVecD;
2050  double* mid = mInvMatD;
2051  mTaperF = static_cast<float>(taper);
2052  mTaperD = taper;
2053  for (int i = 0; i < 3; ++i) {
2054  *vd++ = mat[3][i]; //translation
2055  *vf++ = static_cast<float>(mat[3][i]);
2056  for (int j = 0; j < 3; ++j) {
2057  *md++ = mat[j][i]; //transposed
2058  *mid++ = invMat[j][i];
2059  *mf++ = static_cast<float>(mat[j][i]);
2060  *mif++ = static_cast<float>(invMat[j][i]);
2061  }
2062  }
2063 }
2064 
2065 // ----------------------------> GridBlindMetaData <--------------------------------------
2066 
2067 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) GridBlindMetaData
2068 {
2069  static const int MaxNameSize = 256;// due to NULL termination the maximum length is one less!
2070  int64_t mByteOffset; // byte offset to the blind data, relative to the GridData.
2071  uint64_t mElementCount; // number of elements, e.g. point count
2072  uint32_t mFlags; // flags
2073  GridBlindDataSemantic mSemantic; // semantic meaning of the data.
2075  GridType mDataType; // 4 bytes
2076  char mName[MaxNameSize];// note this include the NULL termination
2077 
2078  /// @brief return memory usage in bytes for the class (note this computes for all blindMetaData structures.)
2079  __hostdev__ static uint64_t memUsage(uint64_t blindDataCount = 0)
2080  {
2081  return blindDataCount * sizeof(GridBlindMetaData);
2082  }
2083 
2084  __hostdev__ void setBlindData(void *ptr) { mByteOffset = PtrDiff(ptr, this); }
2085 
2086  template <typename T>
2087  __hostdev__ const T* getBlindData() const { return PtrAdd<T>(this, mByteOffset); }
2088 
2089 }; // GridBlindMetaData
2090 
2091 // ----------------------------> NodeTrait <--------------------------------------
2092 
2093 /// @brief Struct to derive node type from its level in a given
2094 /// grid, tree or root while perserving constness
2095 template<typename GridOrTreeOrRootT, int LEVEL>
2096 struct NodeTrait;
2097 
2098 // Partial template specialization of above Node struct
2099 template<typename GridOrTreeOrRootT>
2100 struct NodeTrait<GridOrTreeOrRootT, 0>
2101 {
2102  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2103  using Type = typename GridOrTreeOrRootT::LeafNodeType;
2104  using type = typename GridOrTreeOrRootT::LeafNodeType;
2105 };
2106 template<typename GridOrTreeOrRootT>
2107 struct NodeTrait<const GridOrTreeOrRootT, 0>
2108 {
2109  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2110  using Type = const typename GridOrTreeOrRootT::LeafNodeType;
2111  using type = const typename GridOrTreeOrRootT::LeafNodeType;
2112 };
2113 
2114 template<typename GridOrTreeOrRootT>
2115 struct NodeTrait<GridOrTreeOrRootT, 1>
2116 {
2117  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2118  using Type = typename GridOrTreeOrRootT::RootType::ChildNodeType::ChildNodeType;
2119  using type = typename GridOrTreeOrRootT::RootType::ChildNodeType::ChildNodeType;
2120 };
2121 template<typename GridOrTreeOrRootT>
2122 struct NodeTrait<const GridOrTreeOrRootT, 1>
2123 {
2124  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2125  using Type = const typename GridOrTreeOrRootT::RootType::ChildNodeType::ChildNodeType;
2126  using type = const typename GridOrTreeOrRootT::RootType::ChildNodeType::ChildNodeType;
2127 };
2128 template<typename GridOrTreeOrRootT>
2129 struct NodeTrait<GridOrTreeOrRootT, 2>
2130 {
2131  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2132  using Type = typename GridOrTreeOrRootT::RootType::ChildNodeType;
2133  using type = typename GridOrTreeOrRootT::RootType::ChildNodeType;
2134 };
2135 template<typename GridOrTreeOrRootT>
2136 struct NodeTrait<const GridOrTreeOrRootT, 2>
2137 {
2138  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2139  using Type = const typename GridOrTreeOrRootT::RootType::ChildNodeType;
2140  using type = const typename GridOrTreeOrRootT::RootType::ChildNodeType;
2141 };
2142 template<typename GridOrTreeOrRootT>
2143 struct NodeTrait<GridOrTreeOrRootT, 3>
2144 {
2145  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2146  using Type = typename GridOrTreeOrRootT::RootType;
2147  using type = typename GridOrTreeOrRootT::RootType;
2148 };
2149 
2150 template<typename GridOrTreeOrRootT>
2151 struct NodeTrait<const GridOrTreeOrRootT, 3>
2152 {
2153  static_assert(GridOrTreeOrRootT::RootType::LEVEL == 3, "Tree depth is not supported");
2154  using Type = const typename GridOrTreeOrRootT::RootType;
2155  using type = const typename GridOrTreeOrRootT::RootType;
2156 };
2157 
2158 // ----------------------------> Grid <--------------------------------------
2159 
2160 /*
2161  The following class and comment is for internal use only
2162 
2163  Memory layout:
2164 
2165  Grid -> 39 x double (world bbox and affine transformation)
2166  Tree -> Root 3 x ValueType + int32_t + N x Tiles (background,min,max,tileCount + tileCount x Tiles)
2167 
2168  N2 upper InternalNodes each with 2 bit masks, N2 tiles, and min/max values
2169 
2170  N1 lower InternalNodes each with 2 bit masks, N1 tiles, and min/max values
2171 
2172  N0 LeafNodes each with a bit mask, N0 ValueTypes and min/max
2173 
2174  Example layout: ("---" implies it has a custom offset, "..." implies zero or more)
2175  [GridData][TreeData]---[RootData][ROOT TILES...]---[NodeData<5>]---[ModeData<4>]---[LeafData<3>]---[BLINDMETA...]---[BLIND0]---[BLIND1]---etc.
2176 */
2177 
2178 /// @brief Struct with all the member data of the Grid (useful during serialization of an openvdb grid)
2179 ///
2180 /// @note The transform is assumed to be affine (so linear) and have uniform scale! So frustum transforms
2181 /// and non-uniform scaling are not supported (primarily because they complicate ray-tracing in index space)
2182 ///
2183 /// @note No client code should (or can) interface with this struct so it can safely be ignored!
2184 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) GridData
2185 {// sizeof(GridData) = 672B
2186  static const int MaxNameSize = 256;// due to NULL termination the maximum length is one less
2187  uint64_t mMagic; // 8B magic to validate it is valid grid data.
2188  uint64_t mChecksum; // 8B. Checksum of grid buffer.
2189  Version mVersion;// 4B major, minor, and patch version numbers
2190  uint32_t mFlags; // 4B. flags for grid.
2191  uint32_t mGridIndex;// 4B. Index of this grid in the buffer
2192  uint32_t mGridCount; // 4B. Total number of grids in the buffer
2193  uint64_t mGridSize; // 8B. byte count of this entire grid occupied in the buffer.
2194  char mGridName[MaxNameSize]; // 256B
2195  Map mMap; // 264B. affine transformation between index and world space in both single and double precision
2196  BBox<Vec3R> mWorldBBox; // 48B. floating-point AABB of active values in WORLD SPACE (2 x 3 doubles)
2197  Vec3R mVoxelSize; // 24B. size of a voxel in world units
2200  int64_t mBlindMetadataOffset; // 8B. offset of GridBlindMetaData structures that follow this grid.
2201  uint32_t mBlindMetadataCount; // 4B. count of GridBlindMetaData structures that follow this grid.
2202 
2203 
2204  // Set and unset various bit flags
2205  __hostdev__ void setFlagsOff() { mFlags = uint32_t(0); }
2206  __hostdev__ void setMinMaxOn(bool on = true)
2207  {
2208  if (on) {
2209  mFlags |= static_cast<uint32_t>(GridFlags::HasMinMax);
2210  } else {
2211  mFlags &= ~static_cast<uint32_t>(GridFlags::HasMinMax);
2212  }
2213  }
2214  __hostdev__ void setBBoxOn(bool on = true)
2215  {
2216  if (on) {
2217  mFlags |= static_cast<uint32_t>(GridFlags::HasBBox);
2218  } else {
2219  mFlags &= ~static_cast<uint32_t>(GridFlags::HasBBox);
2220  }
2221  }
2222  __hostdev__ void setLongGridNameOn(bool on = true)
2223  {
2224  if (on) {
2225  mFlags |= static_cast<uint32_t>(GridFlags::HasLongGridName);
2226  } else {
2227  mFlags &= ~static_cast<uint32_t>(GridFlags::HasLongGridName);
2228  }
2229  }
2230  __hostdev__ void setAverageOn(bool on = true)
2231  {
2232  if (on) {
2233  mFlags |= static_cast<uint32_t>(GridFlags::HasAverage);
2234  } else {
2235  mFlags &= ~static_cast<uint32_t>(GridFlags::HasAverage);
2236  }
2237  }
2238  __hostdev__ void setStdDeviationOn(bool on = true)
2239  {
2240  if (on) {
2241  mFlags |= static_cast<uint32_t>(GridFlags::HasStdDeviation);
2242  } else {
2243  mFlags &= ~static_cast<uint32_t>(GridFlags::HasStdDeviation);
2244  }
2245  }
2246  __hostdev__ void setBreadthFirstOn(bool on = true)
2247  {
2248  if (on) {
2249  mFlags |= static_cast<uint32_t>(GridFlags::IsBreadthFirst);
2250  } else {
2251  mFlags &= ~static_cast<uint32_t>(GridFlags::IsBreadthFirst);
2252  }
2253  }
2254 
2255  // Affine transformations based on double precision
2256  template<typename Vec3T>
2257  __hostdev__ Vec3T applyMap(const Vec3T& xyz) const { return mMap.applyMap(xyz); } // Pos: index -> world
2258  template<typename Vec3T>
2259  __hostdev__ Vec3T applyInverseMap(const Vec3T& xyz) const { return mMap.applyInverseMap(xyz); } // Pos: world -> index
2260  template<typename Vec3T>
2261  __hostdev__ Vec3T applyJacobian(const Vec3T& xyz) const { return mMap.applyJacobian(xyz); } // Dir: index -> world
2262  template<typename Vec3T>
2263  __hostdev__ Vec3T applyInverseJacobian(const Vec3T& xyz) const { return mMap.applyInverseJacobian(xyz); } // Dir: world -> index
2264  template<typename Vec3T>
2265  __hostdev__ Vec3T applyIJT(const Vec3T& xyz) const { return mMap.applyIJT(xyz); }
2266  // Affine transformations based on single precision
2267  template<typename Vec3T>
2268  __hostdev__ Vec3T applyMapF(const Vec3T& xyz) const { return mMap.applyMapF(xyz); } // Pos: index -> world
2269  template<typename Vec3T>
2270  __hostdev__ Vec3T applyInverseMapF(const Vec3T& xyz) const { return mMap.applyInverseMapF(xyz); } // Pos: world -> index
2271  template<typename Vec3T>
2272  __hostdev__ Vec3T applyJacobianF(const Vec3T& xyz) const { return mMap.applyJacobianF(xyz); } // Dir: index -> world
2273  template<typename Vec3T>
2274  __hostdev__ Vec3T applyInverseJacobianF(const Vec3T& xyz) const { return mMap.applyInverseJacobianF(xyz); } // Dir: world -> index
2275  template<typename Vec3T>
2276  __hostdev__ Vec3T applyIJTF(const Vec3T& xyz) const { return mMap.applyIJTF(xyz); }
2277 
2278  // @brief Return a non-const void pointer to the tree
2279  __hostdev__ void* treePtr() { return this + 1; }
2280 
2281  // @brief Return a const void pointer to the tree
2282  __hostdev__ const void* treePtr() const { return this + 1; }
2283 
2284  /// @brief Returns a const reference to the blindMetaData at the specified linear offset.
2285  ///
2286  /// @warning The linear offset is assumed to be in the valid range
2288  {
2289  NANOVDB_ASSERT(n < mBlindMetadataCount);
2290  return PtrAdd<GridBlindMetaData>(this, mBlindMetadataOffset) + n;
2291  }
2292 
2293 }; // GridData
2294 
2295 // Forward declaration of accelerated random access class
2296 template <typename BuildT, int LEVEL0 = -1, int LEVEL1 = -1, int LEVEL2 = -1>
2298 
2299 template <typename BuildT>
2301 
2302 /// @brief Highest level of the data structure. Contains a tree and a world->index
2303 /// transform (that currently only supports uniform scaling and translation).
2304 ///
2305 /// @note This the API of this class to interface with client code
2306 template<typename TreeT>
2307 class Grid : private GridData
2308 {
2309 public:
2310  using TreeType = TreeT;
2311  using RootType = typename TreeT::RootType;
2313  using ValueType = typename TreeT::ValueType;
2314  using BuildType = typename TreeT::BuildType;// in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
2315  using CoordType = typename TreeT::CoordType;
2317 
2318  /// @brief Disallow constructions, copy and assignment
2319  ///
2320  /// @note Only a Serializer, defined elsewhere, can instantiate this class
2321  Grid(const Grid&) = delete;
2322  Grid& operator=(const Grid&) = delete;
2323  ~Grid() = delete;
2324 
2325  __hostdev__ Version version() const { return DataType::mVersion; }
2326 
2327  __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }
2328 
2329  __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }
2330 
2331  /// @brief Return memory usage in bytes for this class only.
2332  __hostdev__ static uint64_t memUsage() { return sizeof(GridData); }
2333 
2334  /// @brief Return the memory footprint of the entire grid, i.e. including all nodes and blind data
2335  __hostdev__ uint64_t gridSize() const { return DataType::mGridSize; }
2336 
2337  /// @brief Return index of this grid in the buffer
2338  __hostdev__ uint32_t gridIndex() const { return DataType::mGridIndex; }
2339 
2340  /// @brief Return total number of grids in the buffer
2341  __hostdev__ uint32_t gridCount() const { return DataType::mGridCount; }
2342 
2343  /// @brief Return a const reference to the tree
2344  __hostdev__ const TreeT& tree() const { return *reinterpret_cast<const TreeT*>(this->treePtr()); }
2345 
2346  /// @brief Return a non-const reference to the tree
2347  __hostdev__ TreeT& tree() { return *reinterpret_cast<TreeT*>(this->treePtr()); }
2348 
2349  /// @brief Return a new instance of a ReadAccessor used to access values in this grid
2350  __hostdev__ AccessorType getAccessor() const { return AccessorType(this->tree().root()); }
2351 
2352  /// @brief Return a const reference to the size of a voxel in world units
2353  __hostdev__ const Vec3R& voxelSize() const { return DataType::mVoxelSize; }
2354 
2355  /// @brief Return a const reference to the Map for this grid
2356  __hostdev__ const Map& map() const { return DataType::mMap; }
2357 
2358  /// @brief world to index space transformation
2359  template<typename Vec3T>
2360  __hostdev__ Vec3T worldToIndex(const Vec3T& xyz) const { return this->applyInverseMap(xyz); }
2361 
2362  /// @brief index to world space transformation
2363  template<typename Vec3T>
2364  __hostdev__ Vec3T indexToWorld(const Vec3T& xyz) const { return this->applyMap(xyz); }
2365 
2366  /// @brief transformation from index space direction to world space direction
2367  /// @warning assumes dir to be normalized
2368  template<typename Vec3T>
2369  __hostdev__ Vec3T indexToWorldDir(const Vec3T& dir) const { return this->applyJacobian(dir); }
2370 
2371  /// @brief transformation from world space direction to index space direction
2372  /// @warning assumes dir to be normalized
2373  template<typename Vec3T>
2374  __hostdev__ Vec3T worldToIndexDir(const Vec3T& dir) const { return this->applyInverseJacobian(dir); }
2375 
2376  /// @brief transform the gradient from index space to world space.
2377  /// @details Applies the inverse jacobian transform map.
2378  template<typename Vec3T>
2379  __hostdev__ Vec3T indexToWorldGrad(const Vec3T& grad) const { return this->applyIJT(grad); }
2380 
2381  /// @brief world to index space transformation
2382  template<typename Vec3T>
2383  __hostdev__ Vec3T worldToIndexF(const Vec3T& xyz) const { return this->applyInverseMapF(xyz); }
2384 
2385  /// @brief index to world space transformation
2386  template<typename Vec3T>
2387  __hostdev__ Vec3T indexToWorldF(const Vec3T& xyz) const { return this->applyMapF(xyz); }
2388 
2389  /// @brief transformation from index space direction to world space direction
2390  /// @warning assumes dir to be normalized
2391  template<typename Vec3T>
2392  __hostdev__ Vec3T indexToWorldDirF(const Vec3T& dir) const { return this->applyJacobianF(dir); }
2393 
2394  /// @brief transformation from world space direction to index space direction
2395  /// @warning assumes dir to be normalized
2396  template<typename Vec3T>
2397  __hostdev__ Vec3T worldToIndexDirF(const Vec3T& dir) const { return this->applyInverseJacobianF(dir); }
2398 
2399  /// @brief Transforms the gradient from index space to world space.
2400  /// @details Applies the inverse jacobian transform map.
2401  template<typename Vec3T>
2402  __hostdev__ Vec3T indexToWorldGradF(const Vec3T& grad) const { return DataType::applyIJTF(grad); }
2403 
2404  /// @brief Computes a AABB of active values in world space
2405  __hostdev__ const BBox<Vec3R>& worldBBox() const { return DataType::mWorldBBox; }
2406 
2407  /// @brief Computes a AABB of active values in index space
2408  ///
2409  /// @note This method is returning a floating point bounding box and not a CoordBBox. This makes
2410  /// it more useful for clipping rays.
2411  __hostdev__ const BBox<CoordType>& indexBBox() const { return this->tree().bbox(); }
2412 
2413  /// @brief Return the total number of active voxels in this tree.
2414  __hostdev__ uint64_t activeVoxelCount() const { return this->tree().activeVoxelCount(); }
2415 
2416  /// @brief Methods related to the classification of this grid
2417  __hostdev__ bool isValid() const { return DataType::mMagic == NANOVDB_MAGIC_NUMBER; }
2418  __hostdev__ const GridType& gridType() const { return DataType::mGridType; }
2419  __hostdev__ const GridClass& gridClass() const { return DataType::mGridClass; }
2420  __hostdev__ bool isLevelSet() const { return DataType::mGridClass == GridClass::LevelSet; }
2421  __hostdev__ bool isFogVolume() const { return DataType::mGridClass == GridClass::FogVolume; }
2422  __hostdev__ bool isStaggered() const { return DataType::mGridClass == GridClass::Staggered; }
2423  __hostdev__ bool isPointIndex() const { return DataType::mGridClass == GridClass::PointIndex; }
2424  __hostdev__ bool isPointData() const { return DataType::mGridClass == GridClass::PointData; }
2425  __hostdev__ bool isMask() const { return DataType::mGridClass == GridClass::Topology; }
2426  __hostdev__ bool isUnknown() const { return DataType::mGridClass == GridClass::Unknown; }
2427  __hostdev__ bool hasMinMax() const { return DataType::mFlags & static_cast<uint32_t>(GridFlags::HasMinMax); }
2428  __hostdev__ bool hasBBox() const { return DataType::mFlags & static_cast<uint32_t>(GridFlags::HasBBox); }
2429  __hostdev__ bool hasLongGridName() const { return DataType::mFlags & static_cast<uint32_t>(GridFlags::HasLongGridName); }
2430  __hostdev__ bool hasAverage() const { return DataType::mFlags & static_cast<uint32_t>(GridFlags::HasAverage); }
2431  __hostdev__ bool hasStdDeviation() const { return DataType::mFlags & static_cast<uint32_t>(GridFlags::HasStdDeviation); }
2432  __hostdev__ bool isBreadthFirst() const { return DataType::mFlags & static_cast<uint32_t>(GridFlags::IsBreadthFirst); }
2433 
2434  /// @brief return true if the specified node type is layed out breadth-first in memory and has a fixed size.
2435  /// This allows for sequential access to the nodes.
2436  template <typename NodeT>
2437  __hostdev__ bool isSequential() const { return NodeT::FIXED_SIZE && this->isBreadthFirst(); }
2438 
2439  /// @brief return true if the specified node level is layed out breadth-first in memory and has a fixed size.
2440  /// This allows for sequential access to the nodes.
2441  template <int LEVEL>
2442  __hostdev__ bool isSequential() const { return NodeTrait<TreeT,LEVEL>::type::FIXED_SIZE && this->isBreadthFirst(); }
2443 
2444  /// @brief Return a c-string with the name of this grid
2445  __hostdev__ const char* gridName() const
2446  {
2447  if (this->hasLongGridName()) {
2448  const auto &metaData = this->blindMetaData(DataType::mBlindMetadataCount-1);// always the last
2449  NANOVDB_ASSERT(metaData.mDataClass == GridBlindDataClass::GridName);
2450  return metaData.template getBlindData<const char>();
2451  }
2452  return DataType::mGridName;
2453  }
2454 
2455  /// @brief Return a c-string with the name of this grid, truncated to 255 characters
2456  __hostdev__ const char* shortGridName() const { return DataType::mGridName; }
2457 
2458  /// @brief Return checksum of the grid buffer.
2459  __hostdev__ uint64_t checksum() const { return DataType::mChecksum; }
2460 
2461  /// @brief Return true if this grid is empty, i.e. contains no values or nodes.
2462  __hostdev__ bool isEmpty() const { return this->tree().isEmpty(); }
2463 
2464  /// @brief Return the count of blind-data encoded in this grid
2465  __hostdev__ int blindDataCount() const { return DataType::mBlindMetadataCount; }
2466 
2467  /// @brief Return the index of the blind data with specified semantic if found, otherwise -1.
2468  __hostdev__ int findBlindDataForSemantic(GridBlindDataSemantic semantic) const;
2469 
2470  /// @brief Returns a const pointer to the blindData at the specified linear offset.
2471  ///
2472  /// @warning Point might be NULL and the linear offset is assumed to be in the valid range
2473  __hostdev__ const void* blindData(uint32_t n) const
2474  {
2475  if (DataType::mBlindMetadataCount == 0) {
2476  return nullptr;
2477  }
2478  NANOVDB_ASSERT(n < DataType::mBlindMetadataCount);
2479  return this->blindMetaData(n).template getBlindData<void>();
2480  }
2481 
2482  __hostdev__ const GridBlindMetaData& blindMetaData(int n) const { return *DataType::blindMetaData(n); }
2483 
2484 private:
2485  static_assert(sizeof(GridData) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(GridData) is misaligned");
2486 }; // Class Grid
2487 
2488 template<typename TreeT>
2490 {
2491  for (uint32_t i = 0, n = this->blindDataCount(); i < n; ++i)
2492  if (this->blindMetaData(i).mSemantic == semantic)
2493  return int(i);
2494  return -1;
2495 }
2496 
2497 // ----------------------------> Tree <--------------------------------------
2498 
2499 template<int ROOT_LEVEL = 3>
2500 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) TreeData
2501 {// sizeof(TreeData<3>) == 64B
2502  static_assert(ROOT_LEVEL == 3, "Root level is assumed to be three");
2503  uint64_t mNodeOffset[4];//32B, byte offset from this tree to first leaf, lower, upper and root node
2504  uint32_t mNodeCount[3];// 12B, total number of nodes of type: leaf, lower internal, upper internal
2505  uint32_t mTileCount[3];// 12B, total number of tiles of type: leaf, lower internal, upper internal (node, only active tiles!)
2506  uint64_t mVoxelCount;// 8B, total number of active voxels in the root and all its child nodes.
2507 
2508  template <typename RootT>
2509  __hostdev__ void setRoot(const RootT* root) { mNodeOffset[3] = PtrDiff(root, this); }
2510  template <typename RootT>
2511  __hostdev__ RootT* getRoot() { return PtrAdd<RootT>(this, mNodeOffset[3]); }
2512  template <typename RootT>
2513  __hostdev__ const RootT* getRoot() const { return PtrAdd<RootT>(this, mNodeOffset[3]); }
2514 
2515  template <typename NodeT>
2516  __hostdev__ void setFirstNode(const NodeT* node)
2517  {
2518  mNodeOffset[NodeT::LEVEL] = node ? PtrDiff(node, this) : 0;
2519  }
2520 };
2521 
2522 // ----------------------------> GridTree <--------------------------------------
2523 
2524 /// @brief defines a tree type from a grid type while perserving constness
2525 template<typename GridT>
2526 struct GridTree
2527 {
2528  using Type = typename GridT::TreeType;
2529  using type = typename GridT::TreeType;
2530 };
2531 template<typename GridT>
2532 struct GridTree<const GridT>
2533 {
2534  using Type = const typename GridT::TreeType;
2535  using type = const typename GridT::TreeType;
2536 };
2537 
2538 // ----------------------------> Tree <--------------------------------------
2539 
2540 /// @brief VDB Tree, which is a thin wrapper around a RootNode.
2541 template<typename RootT>
2542 class Tree : private TreeData<RootT::LEVEL>
2543 {
2544  static_assert(RootT::LEVEL == 3, "Tree depth is not supported");
2545  static_assert(RootT::ChildNodeType::LOG2DIM == 5, "Tree configuration is not supported");
2546  static_assert(RootT::ChildNodeType::ChildNodeType::LOG2DIM == 4, "Tree configuration is not supported");
2547  static_assert(RootT::LeafNodeType::LOG2DIM == 3, "Tree configuration is not supported");
2548 
2549 public:
2551  using RootType = RootT;
2552  using LeafNodeType = typename RootT::LeafNodeType;
2553  using ValueType = typename RootT::ValueType;
2554  using BuildType = typename RootT::BuildType;// in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
2555  using CoordType = typename RootT::CoordType;
2557 
2558  using Node3 = RootT;
2559  using Node2 = typename RootT::ChildNodeType;
2560  using Node1 = typename Node2::ChildNodeType;
2562 
2563  /// @brief This class cannot be constructed or deleted
2564  Tree() = delete;
2565  Tree(const Tree&) = delete;
2566  Tree& operator=(const Tree&) = delete;
2567  ~Tree() = delete;
2568 
2569  __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }
2570 
2571  __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }
2572 
2573  /// @brief return memory usage in bytes for the class
2574  __hostdev__ static uint64_t memUsage() { return sizeof(DataType); }
2575 
2576  __hostdev__ RootT& root() { return *DataType::template getRoot<RootT>(); }
2577 
2578  __hostdev__ const RootT& root() const { return *DataType::template getRoot<RootT>(); }
2579 
2580  __hostdev__ AccessorType getAccessor() const { return AccessorType(this->root()); }
2581 
2582  /// @brief Return the value of the given voxel (regardless of state or location in the tree.)
2583  __hostdev__ ValueType getValue(const CoordType& ijk) const { return this->root().getValue(ijk); }
2584 
2585  /// @brief Return the active state of the given voxel (regardless of state or location in the tree.)
2586  __hostdev__ bool isActive(const CoordType& ijk) const { return this->root().isActive(ijk); }
2587 
2588  /// @brief Return true if this tree is empty, i.e. contains no values or nodes
2589  __hostdev__ bool isEmpty() const { return this->root().isEmpty(); }
2590 
2591  /// @brief Combines the previous two methods in a single call
2592  __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const { return this->root().probeValue(ijk, v); }
2593 
2594  /// @brief Return a const reference to the background value.
2595  __hostdev__ const ValueType& background() const { return this->root().background(); }
2596 
2597  /// @brief Sets the extrema values of all the active values in this tree, i.e. in all nodes of the tree
2598  __hostdev__ void extrema(ValueType& min, ValueType& max) const;
2599 
2600  /// @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
2601  __hostdev__ const BBox<CoordType>& bbox() const { return this->root().bbox(); }
2602 
2603  /// @brief Return the total number of active voxels in this tree.
2604  __hostdev__ uint64_t activeVoxelCount() const { return DataType::mVoxelCount; }
2605 
2606  /// @brief Return the total number of active tiles at the specified level of the tree.
2607  ///
2608  /// @details n = 0 corresponds to leaf level tiles.
2609  __hostdev__ const uint32_t& activeTileCount(uint32_t n) const
2610  {
2611  NANOVDB_ASSERT(n < 3);
2612  return DataType::mTileCount[n];
2613  }
2614 
2615  template<typename NodeT>
2616  __hostdev__ uint32_t nodeCount() const
2617  {
2618  static_assert(NodeT::LEVEL < 3, "Invalid NodeT");
2619  return DataType::mNodeCount[NodeT::LEVEL];
2620  }
2621 
2622  __hostdev__ uint32_t nodeCount(int level) const
2623  {
2624  NANOVDB_ASSERT(level < 3);
2625  return DataType::mNodeCount[level];
2626  }
2627 
2628  /// @brief return a pointer to the first node of the specified type
2629  ///
2630  /// @warning Note it may return NULL if no nodes exist
2631  template <typename NodeT>
2633  {
2634  const uint64_t offset = DataType::mNodeOffset[NodeT::LEVEL];
2635  return offset>0 ? PtrAdd<NodeT>(this, offset) : nullptr;
2636  }
2637 
2638  /// @brief return a const pointer to the first node of the specified type
2639  ///
2640  /// @warning Note it may return NULL if no nodes exist
2641  template <typename NodeT>
2642  __hostdev__ const NodeT* getFirstNode() const
2643  {
2644  const uint64_t offset = DataType::mNodeOffset[NodeT::LEVEL];
2645  return offset>0 ? PtrAdd<NodeT>(this, offset) : nullptr;
2646  }
2647 
2648  /// @brief return a pointer to the first node at the specified level
2649  ///
2650  /// @warning Note it may return NULL if no nodes exist
2651  template <int LEVEL>
2654  {
2655  return this->template getFirstNode<typename NodeTrait<RootT,LEVEL>::type>();
2656  }
2657 
2658  /// @brief return a const pointer to the first node of the specified level
2659  ///
2660  /// @warning Note it may return NULL if no nodes exist
2661  template <int LEVEL>
2664  {
2665  return this->template getFirstNode<typename NodeTrait<RootT,LEVEL>::type>();
2666  }
2667 
2668 private:
2669  static_assert(sizeof(DataType) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(TreeData) is misaligned");
2670 
2671 }; // Tree class
2672 
2673 template<typename RootT>
2675 {
2676  min = this->root().minimum();
2677  max = this->root().maximum();
2678 }
2679 
2680 // --------------------------> RootNode <------------------------------------
2681 
2682 /// @brief Struct with all the member data of the RootNode (useful during serialization of an openvdb RootNode)
2683 ///
2684 /// @note No client code should (or can) interface with this struct so it can safely be ignored!
2685 template<typename ChildT>
2686 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) RootData
2687 {
2688  using ValueT = typename ChildT::ValueType;
2689  using BuildT = typename ChildT::BuildType;// in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
2690  using CoordT = typename ChildT::CoordType;
2691  using StatsT = typename ChildT::FloatType;
2692  static constexpr bool FIXED_SIZE = false;
2693 
2694  /// @brief Return a key based on the coordinates of a voxel
2695 #ifdef USE_SINGLE_ROOT_KEY
2696  using KeyT = uint64_t;
2697  template <typename CoordType>
2698  __hostdev__ static KeyT CoordToKey(const CoordType& ijk)
2699  {
2700  static_assert(sizeof(CoordT) == sizeof(CoordType), "Mismatching sizeof");
2701  static_assert(32 - ChildT::TOTAL <= 21, "Cannot use 64 bit root keys");
2702  return (KeyT(uint32_t(ijk[2]) >> ChildT::TOTAL)) | // z is the lower 21 bits
2703  (KeyT(uint32_t(ijk[1]) >> ChildT::TOTAL) << 21) | // y is the middle 21 bits
2704  (KeyT(uint32_t(ijk[0]) >> ChildT::TOTAL) << 42); // x is the upper 21 bits
2705  }
2706  __hostdev__ static CoordT KeyToCoord(const KeyT& key)
2707  {
2708  static constexpr uint64_t MASK = (1u << 21) - 1;
2709  return CoordT(((key >> 42) & MASK) << ChildT::TOTAL,
2710  ((key >> 21) & MASK) << ChildT::TOTAL,
2711  (key & MASK) << ChildT::TOTAL);
2712  }
2713 #else
2714  using KeyT = CoordT;
2715  __hostdev__ static KeyT CoordToKey(const CoordT& ijk) { return ijk & ~ChildT::MASK; }
2716  __hostdev__ static CoordT KeyToCoord(const KeyT& key) { return key; }
2717 #endif
2718  BBox<CoordT> mBBox; // 24B. AABB if active values in index space.
2719  uint32_t mTableSize; // 4B. number of tiles and child pointers in the root node
2720 
2721  ValueT mBackground; // background value, i.e. value of any unset voxel
2722  ValueT mMinimum; // typically 4B, minmum of all the active values
2723  ValueT mMaximum; // typically 4B, maximum of all the active values
2724  StatsT mAverage; // typically 4B, average of all the active values in this node and its child nodes
2725  StatsT mStdDevi; // typically 4B, standard deviation of all the active values in this node and its child nodes
2726 
2727  struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) Tile
2728  {
2729  template <typename CoordType>
2730  __hostdev__ void setChild(const CoordType& k, const ChildT *ptr, const RootData *data)
2731  {
2732  key = CoordToKey(k);
2733  child = PtrDiff(ptr, data);
2734  }
2735  template <typename CoordType, typename ValueType>
2736  __hostdev__ void setValue(const CoordType& k, bool s, const ValueType &v)
2737  {
2738  key = CoordToKey(k);
2739  state = s;
2740  value = v;
2741  child = 0;
2742  }
2743  __hostdev__ bool isChild() const { return child; }
2744  __hostdev__ CoordT origin() const { return KeyToCoord(key); }
2745  KeyT key; // USE_SINGLE_ROOT_KEY ? 8B : 12B
2746  int64_t child; // 8B. signed byte offset from this node to the child node. 0 means it is a constant tile, so use value.
2747  uint32_t state; // 4B. state of tile value
2748  ValueT value; // value of tile (i.e. no child node)
2749  }; // Tile
2750 
2751  /// @brief Returns a non-const reference to the tile at the specified linear offset.
2752  ///
2753  /// @warning The linear offset is assumed to be in the valid range
2754  __hostdev__ const Tile* tile(uint32_t n) const
2755  {
2756  NANOVDB_ASSERT(n < mTableSize);
2757  return reinterpret_cast<const Tile*>(this + 1) + n;
2758  }
2759  __hostdev__ Tile* tile(uint32_t n)
2760  {
2761  NANOVDB_ASSERT(n < mTableSize);
2762  return reinterpret_cast<Tile*>(this + 1) + n;
2763  }
2764 
2765  /// @brief Returns a const reference to the child node in the specified tile.
2766  ///
2767  /// @warning A child node is assumed to exist in the specified tile
2768  __hostdev__ ChildT* getChild(const Tile* tile)
2769  {
2770  NANOVDB_ASSERT(tile->child);
2771  return PtrAdd<ChildT>(this, tile->child);
2772  }
2773  __hostdev__ const ChildT* getChild(const Tile* tile) const
2774  {
2775  NANOVDB_ASSERT(tile->child);
2776  return PtrAdd<ChildT>(this, tile->child);
2777  }
2778 
2779  __hostdev__ const ValueT& getMin() const { return mMinimum; }
2780  __hostdev__ const ValueT& getMax() const { return mMaximum; }
2781  __hostdev__ const StatsT& average() const { return mAverage; }
2782  __hostdev__ const StatsT& stdDeviation() const { return mStdDevi; }
2783 
2784  __hostdev__ void setMin(const ValueT& v) { mMinimum = v; }
2785  __hostdev__ void setMax(const ValueT& v) { mMaximum = v; }
2786  __hostdev__ void setAvg(const StatsT& v) { mAverage = v; }
2787  __hostdev__ void setDev(const StatsT& v) { mStdDevi = v; }
2788 
2789  /// @brief This class cannot be constructed or deleted
2790  RootData() = delete;
2791  RootData(const RootData&) = delete;
2792  RootData& operator=(const RootData&) = delete;
2793  ~RootData() = delete;
2794 }; // RootData
2795 
2796 /// @brief Top-most node of the VDB tree structure.
2797 template<typename ChildT>
2798 class RootNode : private RootData<ChildT>
2799 {
2800 public:
2802  using LeafNodeType = typename ChildT::LeafNodeType;
2803  using ChildNodeType = ChildT;
2804  using RootType = RootNode<ChildT>;// this allows RootNode to behave like a Tree
2805 
2806  using ValueType = typename DataType::ValueT;
2807  using FloatType = typename DataType::StatsT;
2808  using BuildType = typename DataType::BuildT;// in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
2809 
2810  using CoordType = typename ChildT::CoordType;
2812  using Tile = typename DataType::Tile;
2813  static constexpr bool FIXED_SIZE = DataType::FIXED_SIZE;
2814 
2815  static constexpr uint32_t LEVEL = 1 + ChildT::LEVEL; // level 0 = leaf
2816 
2817  /// @brief This class cannot be constructed or deleted
2818  RootNode() = delete;
2819  RootNode(const RootNode&) = delete;
2820  RootNode& operator=(const RootNode&) = delete;
2821  ~RootNode() = delete;
2822 
2824 
2825  __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }
2826 
2827  __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }
2828 
2829  /// @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
2830  __hostdev__ const BBox<CoordType>& bbox() const { return DataType::mBBox; }
2831 
2832  /// @brief Return the total number of active voxels in the root and all its child nodes.
2833 
2834  /// @brief Return a const reference to the background value, i.e. the value associated with
2835  /// any coordinate location that has not been set explicitly.
2836  __hostdev__ const ValueType& background() const { return DataType::mBackground; }
2837 
2838  /// @brief Return the number of tiles encoded in this root node
2839  __hostdev__ const uint32_t& tileCount() const { return DataType::mTableSize; }
2840 
2841  /// @brief Return a const reference to the minimum active value encoded in this root node and any of its child nodes
2842  __hostdev__ const ValueType& minimum() const { return this->getMin(); }
2843 
2844  /// @brief Return a const reference to the maximum active value encoded in this root node and any of its child nodes
2845  __hostdev__ const ValueType& maximum() const { return this->getMax(); }
2846 
2847  /// @brief Return a const reference to the average of all the active values encoded in this root node and any of its child nodes
2848  __hostdev__ const FloatType& average() const { return DataType::mAverage; }
2849 
2850  /// @brief Return the variance of all the active values encoded in this root node and any of its child nodes
2851  __hostdev__ FloatType variance() const { return DataType::mStdDevi * DataType::mStdDevi; }
2852 
2853  /// @brief Return a const reference to the standard deviation of all the active values encoded in this root node and any of its child nodes
2854  __hostdev__ const FloatType& stdDeviation() const { return DataType::mStdDevi; }
2855 
2856  /// @brief Return the expected memory footprint in bytes with the specified number of tiles
2857  __hostdev__ static uint64_t memUsage(uint32_t tableSize) { return sizeof(RootNode) + tableSize * sizeof(Tile); }
2858 
2859  /// @brief Return the actual memory footprint of this root node
2860  __hostdev__ uint64_t memUsage() const { return sizeof(RootNode) + DataType::mTableSize * sizeof(Tile); }
2861 
2862  /// @brief Return the value of the given voxel
2864  {
2865  if (const Tile* tile = this->findTile(ijk)) {
2866  return tile->isChild() ? this->getChild(tile)->getValue(ijk) : tile->value;
2867  }
2868  return DataType::mBackground;
2869  }
2870 
2871  __hostdev__ bool isActive(const CoordType& ijk) const
2872  {
2873  if (const Tile* tile = this->findTile(ijk)) {
2874  return tile->isChild() ? this->getChild(tile)->isActive(ijk) : tile->state;
2875  }
2876  return false;
2877  }
2878 
2879  /// @brief Return true if this RootNode is empty, i.e. contains no values or nodes
2880  __hostdev__ bool isEmpty() const { return DataType::mTableSize == uint32_t(0); }
2881 
2882  __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const
2883  {
2884  if (const Tile* tile = this->findTile(ijk)) {
2885  if (tile->isChild()) {
2886  const auto *child = this->getChild(tile);
2887  return child->probeValue(ijk, v);
2888  }
2889  v = tile->value;
2890  return tile->state;
2891  }
2892  v = DataType::mBackground;
2893  return false;
2894  }
2895 
2896  __hostdev__ const LeafNodeType* probeLeaf(const CoordType& ijk) const
2897  {
2898  const Tile* tile = this->findTile(ijk);
2899  if (tile && tile->isChild()) {
2900  const auto *child = this->getChild(tile);
2901  return child->probeLeaf(ijk);
2902  }
2903  return nullptr;
2904  }
2905 
2906 private:
2907  static_assert(sizeof(DataType) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(RootData) is misaligned");
2908  static_assert(sizeof(typename DataType::Tile) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(RootData::Tile) is misaligned");
2909 
2910  template<typename, int, int, int>
2911  friend class ReadAccessor;
2912 
2913  template<typename>
2914  friend class Tree;
2915 
2916  /// @brief Private method to find a Tile of this root node by means of binary-search. This is obviously
2917  /// much slower then direct lookup into a linear array (as in the other nodes) which is exactly
2918  /// why it is important to use the ReadAccessor which amortizes this overhead by node caching and
2919  /// inverse tree traversal!
2920  __hostdev__ const Tile* findTile(const CoordType& ijk) const
2921  {
2922  const Tile* tiles = reinterpret_cast<const Tile*>(this + 1);
2923  const auto key = DataType::CoordToKey(ijk);
2924 #if 1 // switch between linear and binary seach
2925  for (uint32_t i = 0; i < DataType::mTableSize; ++i) {
2926  if (tiles[i].key == key) return &tiles[i];
2927  }
2928 #else// do not enable binary search if tiles are not guaranteed to be sorted!!!!!!
2929  // binary-search of pre-sorted elements
2930  int32_t low = 0, high = DataType::mTableSize; // low is inclusive and high is exclusive
2931  while (low != high) {
2932  int mid = low + ((high - low) >> 1);
2933  const Tile* tile = &tiles[mid];
2934  if (tile->key == key) {
2935  return tile;
2936  } else if (tile->key < key) {
2937  low = mid + 1;
2938  } else {
2939  high = mid;
2940  }
2941  }
2942 #endif
2943  return nullptr;
2944  }
2945 
2946  /// @brief Private method to return node information and update a ReadAccessor
2947  template<typename AccT>
2948  __hostdev__ typename AccT::NodeInfo getNodeInfoAndCache(const CoordType& ijk, const AccT& acc) const
2949  {
2950  using NodeInfoT = typename AccT::NodeInfo;
2951  if (const Tile* tile = this->findTile(ijk)) {
2952  if (tile->isChild()) {
2953  const auto *child = this->getChild(tile);
2954  acc.insert(ijk, child);
2955  return child->getNodeInfoAndCache(ijk, acc);
2956  }
2957  return NodeInfoT{LEVEL, ChildT::dim(), tile->value, tile->value, tile->value,
2958  0, tile->origin(), tile->origin() + CoordType(ChildT::DIM)};
2959  }
2960  return NodeInfoT{LEVEL, ChildT::dim(), this->minimum(), this->maximum(),
2961  this->average(), this->stdDeviation(), this->bbox()[0], this->bbox()[1]};
2962  }
2963 
2964  /// @brief Private method to return a voxel value and update a ReadAccessor
2965  template<typename AccT>
2966  __hostdev__ ValueType getValueAndCache(const CoordType& ijk, const AccT& acc) const
2967  {
2968  if (const Tile* tile = this->findTile(ijk)) {
2969  if (tile->isChild()) {
2970  const auto *child = this->getChild(tile);
2971  acc.insert(ijk, child);
2972  return child->getValueAndCache(ijk, acc);
2973  }
2974  return tile->value;
2975  }
2976  return DataType::mBackground;
2977  }
2978 
2979  template<typename AccT>
2980  __hostdev__ bool isActiveAndCache(const CoordType& ijk, const AccT& acc) const
2981  {
2982  const Tile* tile = this->findTile(ijk);
2983  if (tile && tile->isChild()) {
2984  const auto *child = this->getChild(tile);
2985  acc.insert(ijk, child);
2986  return child->isActiveAndCache(ijk, acc);
2987  }
2988  return false;
2989  }
2990 
2991  template<typename AccT>
2992  __hostdev__ bool probeValueAndCache(const CoordType& ijk, ValueType& v, const AccT& acc) const
2993  {
2994  if (const Tile* tile = this->findTile(ijk)) {
2995  if (tile->isChild()) {
2996  const auto *child = this->getChild(tile);
2997  acc.insert(ijk, child);
2998  return child->probeValueAndCache(ijk, v, acc);
2999  }
3000  v = tile->value;
3001  return tile->state;
3002  }
3003  v = DataType::mBackground;
3004  return false;
3005  }
3006 
3007  template<typename AccT>
3008  __hostdev__ const LeafNodeType* probeLeafAndCache(const CoordType& ijk, const AccT& acc) const
3009  {
3010  const Tile* tile = this->findTile(ijk);
3011  if (tile && tile->isChild()) {
3012  const auto *child = this->getChild(tile);
3013  acc.insert(ijk, child);
3014  return child->probeLeafAndCache(ijk, acc);
3015  }
3016  return nullptr;
3017  }
3018 
3019  template<typename RayT, typename AccT>
3020  __hostdev__ uint32_t getDimAndCache(const CoordType& ijk, const RayT& ray, const AccT& acc) const
3021  {
3022  if (const Tile* tile = this->findTile(ijk)) {
3023  if (tile->isChild()) {
3024  const auto *child = this->getChild(tile);
3025  acc.insert(ijk, child);
3026  return child->getDimAndCache(ijk, ray, acc);
3027  }
3028  return 1 << ChildT::TOTAL; //tile value
3029  }
3030  return ChildNodeType::dim(); // background
3031  }
3032 }; // RootNode class
3033 
3034 // After the RootNode the memory layout is assumed to be the sorted Tiles
3035 
3036 // --------------------------> InternalNode <------------------------------------
3037 
3038 /// @brief Struct with all the member data of the InternalNode (useful during serialization of an openvdb InternalNode)
3039 ///
3040 /// @note No client code should (or can) interface with this struct so it can safely be ignored!
3041 template<typename ChildT, uint32_t LOG2DIM>
3042 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) InternalData
3043 {
3044  using ValueT = typename ChildT::ValueType;
3045  using BuildT = typename ChildT::BuildType;// in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
3046  using StatsT = typename ChildT::FloatType;
3047  using CoordT = typename ChildT::CoordType;
3048  using MaskT = typename ChildT::template MaskType<LOG2DIM>;
3049  static constexpr bool FIXED_SIZE = true;
3050 
3051  union Tile
3052  {
3054  int64_t child;//signed 64 bit byte offset relative to the InternalData!!
3055  /// @brief This class cannot be constructed or deleted
3056  Tile() = delete;
3057  Tile(const Tile&) = delete;
3058  Tile& operator=(const Tile&) = delete;
3059  ~Tile() = delete;
3060  };
3061 
3062  BBox<CoordT> mBBox; // 24B. node bounding box. |
3063  uint64_t mFlags; // 8B. node flags. | 32B aligned
3064  MaskT mValueMask; // LOG2DIM(5): 4096B, LOG2DIM(4): 512B | 32B aligned
3065  MaskT mChildMask; // LOG2DIM(5): 4096B, LOG2DIM(4): 512B | 32B aligned
3066 
3067  ValueT mMinimum; // typically 4B
3068  ValueT mMaximum; // typically 4B
3069  StatsT mAverage; // typically 4B, average of all the active values in this node and its child nodes
3070  StatsT mStdDevi; // typically 4B, standard deviation of all the active values in this node and its child nodes
3071  alignas(32) Tile mTable[1u << (3 * LOG2DIM)]; // sizeof(ValueT) x (16*16*16 or 32*32*32)
3072 
3073  __hostdev__ void setChild(uint32_t n, const void *ptr)
3074  {
3075  NANOVDB_ASSERT(mChildMask.isOn(n));
3076  mTable[n].child = PtrDiff(ptr, this);
3077  }
3078 
3079  template <typename ValueT>
3080  __hostdev__ void setValue(uint32_t n, const ValueT &v)
3081  {
3082  NANOVDB_ASSERT(!mChildMask.isOn(n));
3083  mTable[n].value = v;
3084  }
3085 
3086  /// @brief Returns a pointer to the child node at the specifed linear offset.
3087  __hostdev__ ChildT* getChild(uint32_t n)
3088  {
3089  NANOVDB_ASSERT(mChildMask.isOn(n));
3090  return PtrAdd<ChildT>(this, mTable[n].child);
3091  }
3092  __hostdev__ const ChildT* getChild(uint32_t n) const
3093  {
3094  NANOVDB_ASSERT(mChildMask.isOn(n));
3095  return PtrAdd<ChildT>(this, mTable[n].child);
3096  }
3097 
3098  template <typename T>
3099  __hostdev__ void setOrigin(const T& ijk) { mBBox[0] = ijk; }
3100 
3101  __hostdev__ const ValueT& getMin() const { return mMinimum; }
3102  __hostdev__ const ValueT& getMax() const { return mMaximum; }
3103  __hostdev__ const StatsT& average() const { return mAverage; }
3104  __hostdev__ const StatsT& stdDeviation() const { return mStdDevi; }
3105 
3106  __hostdev__ void setMin(const ValueT& v) { mMinimum = v; }
3107  __hostdev__ void setMax(const ValueT& v) { mMaximum = v; }
3108  __hostdev__ void setAvg(const StatsT& v) { mAverage = v; }
3109  __hostdev__ void setDev(const StatsT& v) { mStdDevi = v; }
3110 
3111  /// @brief This class cannot be constructed or deleted
3112  InternalData() = delete;
3113  InternalData(const InternalData&) = delete;
3114  InternalData& operator=(const InternalData&) = delete;
3115  ~InternalData() = delete;
3116 }; // InternalData
3117 
3118 /// @brief Internal nodes of a VDB treedim(),
3119 template<typename ChildT, uint32_t Log2Dim = ChildT::LOG2DIM + 1>
3120 class InternalNode : private InternalData<ChildT, Log2Dim>
3121 {
3122 public:
3124  using ValueType = typename DataType::ValueT;
3125  using FloatType = typename DataType::StatsT;
3126  using BuildType = typename DataType::BuildT; // in rare cases BuildType != ValueType, e.g. then BuildType = ValueMask and ValueType = bool
3127  using LeafNodeType = typename ChildT::LeafNodeType;
3128  using ChildNodeType = ChildT;
3129  using CoordType = typename ChildT::CoordType;
3130  static constexpr bool FIXED_SIZE = DataType::FIXED_SIZE;
3131  template<uint32_t LOG2>
3132  using MaskType = typename ChildT::template MaskType<LOG2>;
3133 
3134  static constexpr uint32_t LOG2DIM = Log2Dim;
3135  static constexpr uint32_t TOTAL = LOG2DIM + ChildT::TOTAL; // dimension in index space
3136  static constexpr uint32_t DIM = 1u << TOTAL; // number of voxels along each axis of this node
3137  static constexpr uint32_t SIZE = 1u << (3 * LOG2DIM); // number of tile values (or child pointers)
3138  static constexpr uint32_t MASK = (1u << TOTAL) - 1u;
3139  static constexpr uint32_t LEVEL = 1 + ChildT::LEVEL; // level 0 = leaf
3140  static constexpr uint64_t NUM_VALUES = uint64_t(1) << (3 * TOTAL); // total voxel count represented by this node
3141 
3142  /// @brief This class cannot be constructed or deleted
3143  InternalNode() = delete;
3144  InternalNode(const InternalNode&) = delete;
3145  InternalNode& operator=(const InternalNode&) = delete;
3146  ~InternalNode() = delete;
3147 
3148  __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }
3149 
3150  __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }
3151 
3152  /// @brief Return the dimension, in voxel units, of this internal node (typically 8*16 or 8*16*32)
3153  __hostdev__ static uint32_t dim() { return 1u << TOTAL; }
3154 
3155  /// @brief Return memory usage in bytes for the class
3156  __hostdev__ static size_t memUsage() { return sizeof(DataType); }
3157 
3158  /// @brief Return a const reference to the bit mask of active voxels in this internal node
3159  __hostdev__ const MaskType<LOG2DIM>& valueMask() const { return DataType::mValueMask; }
3160 
3161  /// @brief Return a const reference to the bit mask of child nodes in this internal node
3162  __hostdev__ const MaskType<LOG2DIM>& childMask() const { return DataType::mChildMask; }
3163 
3164  /// @brief Return the origin in index space of this leaf node
3165  __hostdev__ CoordType origin() const { return DataType::mBBox.min() & ~MASK; }
3166 
3167  /// @brief Return a const reference to the minimum active value encoded in this internal node and any of its child nodes
3168  __hostdev__ const ValueType& minimum() const { return this->getMin(); }
3169 
3170  /// @brief Return a const reference to the maximum active value encoded in this internal node and any of its child nodes
3171  __hostdev__ const ValueType& maximum() const { return this->getMax(); }
3172 
3173  /// @brief Return a const reference to the average of all the active values encoded in this internal node and any of its child nodes
3174  __hostdev__ const FloatType& average() const { return DataType::mAverage; }
3175 
3176  /// @brief Return the variance of all the active values encoded in this internal node and any of its child nodes
3177  __hostdev__ FloatType variance() const { return DataType::mStdDevi*DataType::mStdDevi; }
3178 
3179  /// @brief Return a const reference to the standard deviation of all the active values encoded in this internal node and any of its child nodes
3180  __hostdev__ const FloatType& stdDeviation() const { return DataType::mStdDevi; }
3181 
3182  /// @brief Return a const reference to the bounding box in index space of active values in this internal node and any of its child nodes
3183  __hostdev__ const BBox<CoordType>& bbox() const { return DataType::mBBox; }
3184 
3185  /// @brief Return the value of the given voxel
3187  {
3188  const uint32_t n = CoordToOffset(ijk);
3189  return DataType::mChildMask.isOn(n) ? this->getChild(n)->getValue(ijk) : DataType::mTable[n].value;
3190  }
3191 
3192  __hostdev__ bool isActive(const CoordType& ijk) const
3193  {
3194  const uint32_t n = CoordToOffset(ijk);
3195  return DataType::mChildMask.isOn(n) ? this->getChild(n)->isActive(ijk) : DataType::mValueMask.isOn(n);
3196  }
3197 
3198  __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const
3199  {
3200  const uint32_t n = CoordToOffset(ijk);
3201  if (DataType::mChildMask.isOn(n))
3202  return this->getChild(n)->probeValue(ijk, v);
3203  v = DataType::mTable[n].value;
3204  return DataType::mValueMask.isOn(n);
3205  }
3206 
3207  __hostdev__ const LeafNodeType* probeLeaf(const CoordType& ijk) const
3208  {
3209  const uint32_t n = CoordToOffset(ijk);
3210  if (DataType::mChildMask.isOn(n))
3211  return this->getChild(n)->probeLeaf(ijk);
3212  return nullptr;
3213  }
3214 
3215  /// @brief Return the linear offset corresponding to the given coordinate
3216  __hostdev__ static uint32_t CoordToOffset(const CoordType& ijk)
3217  {
3218 #if 0
3219  return (((ijk[0] & MASK) >> ChildT::TOTAL) << (2 * LOG2DIM)) +
3220  (((ijk[1] & MASK) >> ChildT::TOTAL) << (LOG2DIM)) +
3221  ((ijk[2] & MASK) >> ChildT::TOTAL);
3222 #else
3223  return (((ijk[0] & MASK) >> ChildT::TOTAL) << (2 * LOG2DIM)) |
3224  (((ijk[1] & MASK) >> ChildT::TOTAL) << (LOG2DIM)) |
3225  ((ijk[2] & MASK) >> ChildT::TOTAL);
3226 #endif
3227  }
3228 
3229  /// @return the local coordinate of the n'th tile or child node
3231  {
3232  NANOVDB_ASSERT(n < SIZE);
3233  const uint32_t m = n & ((1 << 2 * LOG2DIM) - 1);
3234  return Coord(n >> 2 * LOG2DIM, m >> LOG2DIM, m & ((1 << LOG2DIM) - 1));
3235  }
3236 
3237  /// @brief modifies local coordinates to global coordinates of a tile or child node
3239  {
3240  ijk <<= ChildT::TOTAL;
3241  ijk += this->origin();
3242  }
3243 
3245  {
3247  this->localToGlobalCoord(ijk);
3248  return ijk;
3249  }
3250 
3251  /// @brief Retrun true if this node or any of its child nodes contain active values
3252  __hostdev__ bool isActive() const
3253  {
3254  return DataType::mFlags & uint32_t(2);
3255  }
3256 
3257 private:
3258  static_assert(sizeof(DataType) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(InternalData) is misaligned");
3259  //static_assert(offsetof(DataType, mTable) % 32 == 0, "InternalData::mTable is misaligned");
3260 
3261  template<typename, int, int, int>
3262  friend class ReadAccessor;
3263 
3264  template<typename>
3265  friend class RootNode;
3266  template<typename, uint32_t>
3267  friend class InternalNode;
3268 
3269  /// @brief Private read access method used by the ReadAccessor
3270  template<typename AccT>
3271  __hostdev__ ValueType getValueAndCache(const CoordType& ijk, const AccT& acc) const
3272  {
3273  const uint32_t n = CoordToOffset(ijk);
3274  if (!DataType::mChildMask.isOn(n))
3275  return DataType::mTable[n].value;
3276  const ChildT* child = this->getChild(n);
3277  acc.insert(ijk, child);
3278  return child->getValueAndCache(ijk, acc);
3279  }
3280 
3281  template<typename AccT>
3282  __hostdev__ typename AccT::NodeInfo getNodeInfoAndCache(const CoordType& ijk, const AccT& acc) const
3283  {
3284  using NodeInfoT = typename AccT::NodeInfo;
3285  const uint32_t n = CoordToOffset(ijk);
3286  if (!DataType::mChildMask.isOn(n)) {
3287  return NodeInfoT{LEVEL, this->dim(), this->minimum(), this->maximum(), this->average(),
3288  this->stdDeviation(), this->bbox()[0], this->bbox()[1]};
3289  }
3290  const ChildT* child = this->getChild(n);
3291  acc.insert(ijk, child);
3292  return child->getNodeInfoAndCache(ijk, acc);
3293  }
3294 
3295  template<typename AccT>
3296  __hostdev__ bool isActiveAndCache(const CoordType& ijk, const AccT& acc) const
3297  {
3298  const uint32_t n = CoordToOffset(ijk);
3299  if (!DataType::mChildMask.isOn(n))
3300  return DataType::mValueMask.isOn(n);
3301  const ChildT* child = this->getChild(n);
3302  acc.insert(ijk, child);
3303  return child->isActiveAndCache(ijk, acc);
3304  }
3305 
3306  template<typename AccT>
3307  __hostdev__ bool probeValueAndCache(const CoordType& ijk, ValueType& v, const AccT& acc) const
3308  {
3309  const uint32_t n = CoordToOffset(ijk);
3310  if (!DataType::mChildMask.isOn(n)) {
3311  v = DataType::mTable[n].value;
3312  return DataType::mValueMask.isOn(n);
3313  }
3314  const ChildT* child = this->getChild(n);
3315  acc.insert(ijk, child);
3316  return child->probeValueAndCache(ijk, v, acc);
3317  }
3318 
3319  template<typename AccT>
3320  __hostdev__ const LeafNodeType* probeLeafAndCache(const CoordType& ijk, const AccT& acc) const
3321  {
3322  const uint32_t n = CoordToOffset(ijk);
3323  if (!DataType::mChildMask.isOn(n))
3324  return nullptr;
3325  const ChildT* child = this->getChild(n);
3326  acc.insert(ijk, child);
3327  return child->probeLeafAndCache(ijk, acc);
3328  }
3329 
3330  template<typename RayT, typename AccT>
3331  __hostdev__ uint32_t getDimAndCache(const CoordType& ijk, const RayT& ray, const AccT& acc) const
3332  {
3333  if (DataType::mFlags & uint32_t(1))
3334  this->dim(); //ship this node if first bit is set
3335  //if (!ray.intersects( this->bbox() )) return 1<<TOTAL;
3336 
3337  const uint32_t n = CoordToOffset(ijk);
3338  if (DataType::mChildMask.isOn(n)) {
3339  const ChildT* child = this->getChild(n);
3340  acc.insert(ijk, child);
3341  return child->getDimAndCache(ijk, ray, acc);
3342  }
3343  return ChildNodeType::dim(); // tile value
3344  }
3345 
3346 }; // InternalNode class
3347 
3348 // --------------------------> LeafNode <------------------------------------
3349 
3350 /// @brief Stuct with all the member data of the LeafNode (useful during serialization of an openvdb LeafNode)
3351 ///
3352 /// @note No client code should (or can) interface with this struct so it can safely be ignored!
3353 template<typename ValueT, typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3354 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData
3355 {
3356  static_assert(sizeof(CoordT) == sizeof(Coord), "Mismatching sizeof");
3357  static_assert(sizeof(MaskT<LOG2DIM>) == sizeof(Mask<LOG2DIM>), "Mismatching sizeof");
3358  using ValueType = ValueT;
3359  using BuildType = ValueT;
3361  using ArrayType = ValueT;// type used for the internal mValue array
3362  static constexpr bool FIXED_SIZE = true;
3363 
3364  CoordT mBBoxMin; // 12B.
3365  uint8_t mBBoxDif[3]; // 3B.
3366  uint8_t mFlags; // 1B.
3367  MaskT<LOG2DIM> mValueMask; // LOG2DIM(3): 64B.
3368 
3369  ValueType mMinimum; // typically 4B
3370  ValueType mMaximum; // typically 4B
3371  FloatType mAverage; // typically 4B, average of all the active values in this node and its child nodes
3372  FloatType mStdDevi; // typically 4B, standard deviation of all the active values in this node and its child nodes
3373  alignas(32) ValueType mValues[1u << 3 * LOG2DIM];
3374 
3375  //__hostdev__ const ValueType* values() const { return mValues; }
3376  __hostdev__ ValueType getValue(uint32_t i) const { return mValues[i]; }
3377  __hostdev__ void setValueOnly(uint32_t offset, const ValueType& value) { mValues[offset] = value; }
3378  __hostdev__ void setValue(uint32_t offset, const ValueType& value)
3379  {
3380  mValueMask.setOn(offset);
3381  mValues[offset] = value;
3382  }
3383 
3384  __hostdev__ ValueType getMin() const { return mMinimum; }
3385  __hostdev__ ValueType getMax() const { return mMaximum; }
3386  __hostdev__ FloatType getAvg() const { return mAverage; }
3387  __hostdev__ FloatType getDev() const { return mStdDevi; }
3388 
3389  __hostdev__ void setMin(const ValueType& v) { mMinimum = v; }
3390  __hostdev__ void setMax(const ValueType& v) { mMaximum = v; }
3391  __hostdev__ void setAvg(const FloatType& v) { mAverage = v; }
3392  __hostdev__ void setDev(const FloatType& v) { mStdDevi = v; }
3393 
3394  template <typename T>
3395  __hostdev__ void setOrigin(const T& ijk) { mBBoxMin = ijk; }
3396 
3397  /// @brief This class cannot be constructed or deleted
3398  LeafData() = delete;
3399  LeafData(const LeafData&) = delete;
3400  LeafData& operator=(const LeafData&) = delete;
3401  ~LeafData() = delete;
3402 }; // LeafData<ValueT>
3403 
3404 /// @brief Base-class for quantized float leaf nodes
3405 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3406 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafFnBase
3407 {
3408  static_assert(sizeof(CoordT) == sizeof(Coord), "Mismatching sizeof");
3409  static_assert(sizeof(MaskT<LOG2DIM>) == sizeof(Mask<LOG2DIM>), "Mismatching sizeof");
3410  using ValueType = float;
3411  using FloatType = float;
3412 
3413  CoordT mBBoxMin; // 12B.
3414  uint8_t mBBoxDif[3]; // 3B.
3415  uint8_t mFlags; // 1B.
3416  MaskT<LOG2DIM> mValueMask; // LOG2DIM(3): 64B.
3417 
3418  float mMinimum; // 4B - minimum of ALL values in this node
3419  float mQuantum; // = (max - min)/15 4B
3420  uint16_t mMin, mMax, mAvg, mDev;// quantized representations of statistics of active values
3421 
3422  void init(float min, float max, uint8_t bitWidth)
3423  {
3424  mMinimum = min;
3425  mQuantum = (max - min)/float((1 << bitWidth)-1);
3426  }
3427 
3428  /// @brief return the quantized minimum of the active values in this node
3429  __hostdev__ float getMin() const { return mMin*mQuantum + mMinimum; }
3430 
3431  /// @brief return the quantized maximum of the active values in this node
3432  __hostdev__ float getMax() const { return mMax*mQuantum + mMinimum; }
3433 
3434  /// @brief return the quantized average of the active values in this node
3435  __hostdev__ float getAvg() const { return mAvg*mQuantum + mMinimum; }
3436  /// @brief return the quantized standard deviation of the active values in this node
3437 
3438  /// @note 0 <= StdDev <= max-min or 0 <= StdDev/(max-min) <= 1
3439  __hostdev__ float getDev() const { return mDev*mQuantum; }
3440 
3441  /// @note min <= X <= max or 0 <= (X-min)/(min-max) <= 1
3442  __hostdev__ void setMin(float min) { mMin = uint16_t((min - mMinimum)/mQuantum + 0.5f); }
3443 
3444  /// @note min <= X <= max or 0 <= (X-min)/(min-max) <= 1
3445  __hostdev__ void setMax(float max) { mMax = uint16_t((max - mMinimum)/mQuantum + 0.5f); }
3446 
3447  /// @note min <= avg <= max or 0 <= (avg-min)/(min-max) <= 1
3448  __hostdev__ void setAvg(float avg) { mAvg = uint16_t((avg - mMinimum)/mQuantum + 0.5f); }
3449 
3450  /// @note 0 <= StdDev <= max-min or 0 <= StdDev/(max-min) <= 1
3451  __hostdev__ void setDev(float dev) { mDev = uint16_t(dev/mQuantum + 0.5f); }
3452 
3453  template <typename T>
3454  __hostdev__ void setOrigin(const T& ijk) { mBBoxMin = ijk; }
3455 };// LeafFnBase
3456 
3457 /// @brief Stuct with all the member data of the LeafNode (useful during serialization of an openvdb LeafNode)
3458 ///
3459 /// @note No client code should (or can) interface with this struct so it can safely be ignored!
3460 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3461 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData<Fp4, CoordT, MaskT, LOG2DIM>
3462  : public LeafFnBase<CoordT, MaskT, LOG2DIM>
3463 {
3465  using BuildType = Fp4;
3466  using ArrayType = uint8_t;// type used for the internal mValue array
3467  static constexpr bool FIXED_SIZE = true;
3468  alignas(32) uint8_t mCode[1u << (3 * LOG2DIM - 1)];
3469 
3470  __hostdev__ static constexpr uint8_t bitWidth() { return 4u; }
3471  __hostdev__ float getValue(uint32_t i) const
3472  {
3473 #if 0
3474  const uint8_t c = mCode[i>>1];
3475  return ( (i&1) ? c >> 4 : c & uint8_t(15) )*BaseT::mQuantum + BaseT::mMinimum;
3476 #else
3477  return ((mCode[i>>1] >> ((i&1)<<2)) & uint8_t(15))*BaseT::mQuantum + BaseT::mMinimum;
3478 #endif
3479  }
3480 
3481  /// @brief This class cannot be constructed or deleted
3482  LeafData() = delete;
3483  LeafData(const LeafData&) = delete;
3484  LeafData& operator=(const LeafData&) = delete;
3485  ~LeafData() = delete;
3486 }; // LeafData<Fp4>
3487 
3488 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3489 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData<Fp8, CoordT, MaskT, LOG2DIM>
3490  : public LeafFnBase<CoordT, MaskT, LOG2DIM>
3491 {
3493  using BuildType = Fp8;
3494  using ArrayType = uint8_t;// type used for the internal mValue array
3495  static constexpr bool FIXED_SIZE = true;
3496  alignas(32) uint8_t mCode[1u << 3 * LOG2DIM];
3497 
3498  __hostdev__ static constexpr uint8_t bitWidth() { return 8u; }
3499  __hostdev__ float getValue(uint32_t i) const
3500  {
3501  return mCode[i]*BaseT::mQuantum + BaseT::mMinimum;// code * (max-min)/255 + min
3502  }
3503  /// @brief This class cannot be constructed or deleted
3504  LeafData() = delete;
3505  LeafData(const LeafData&) = delete;
3506  LeafData& operator=(const LeafData&) = delete;
3507  ~LeafData() = delete;
3508 }; // LeafData<Fp8>
3509 
3510 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3511 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData<Fp16, CoordT, MaskT, LOG2DIM>
3512  : public LeafFnBase<CoordT, MaskT, LOG2DIM>
3513 {
3515  using BuildType = Fp16;
3516  using ArrayType = uint16_t;// type used for the internal mValue array
3517  static constexpr bool FIXED_SIZE = true;
3518  alignas(32) uint16_t mCode[1u << 3 * LOG2DIM];
3519 
3520  __hostdev__ static constexpr uint8_t bitWidth() { return 16u; }
3521  __hostdev__ float getValue(uint32_t i) const
3522  {
3523  return mCode[i]*BaseT::mQuantum + BaseT::mMinimum;// code * (max-min)/65535 + min
3524  }
3525 
3526  /// @brief This class cannot be constructed or deleted
3527  LeafData() = delete;
3528  LeafData(const LeafData&) = delete;
3529  LeafData& operator=(const LeafData&) = delete;
3530  ~LeafData() = delete;
3531 }; // LeafData<Fp16>
3532 
3533 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3534 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData<FpN, CoordT, MaskT, LOG2DIM>
3535  : public LeafFnBase<CoordT, MaskT, LOG2DIM>
3536 {
3538  using BuildType = FpN;
3539  static constexpr bool FIXED_SIZE = false;
3540 
3541  __hostdev__ uint8_t bitWidth() const { return 1 << (BaseT::mFlags >> 5); }// 4,8,16,32 = 2^(2,3,4,5)
3542  __hostdev__ size_t memUsage() const { return sizeof(*this) + this->bitWidth()*64; }
3543  __hostdev__ static size_t memUsage(uint32_t bitWidth) { return 96u + bitWidth*64; }
3544  __hostdev__ float getValue(uint32_t i) const
3545  {
3546 #ifdef NANOVDB_FPN_BRANCHLESS// faster
3547  const int b = BaseT::mFlags >> 5;// b = 0, 1, 2, 3, 4 corresponding to 1, 2, 4, 8, 16 bits
3548 #if 0// use LUT
3549  uint16_t code = reinterpret_cast<const uint16_t*>(this + 1)[i >> (4 - b)];
3550  const static uint8_t shift[5] = {15, 7, 3, 1, 0};
3551  const static uint16_t mask[5] = {1, 3, 15, 255, 65535};
3552  code >>= (i & shift[b]) << b;
3553  code &= mask[b];
3554 #else// no LUT
3555  uint32_t code = reinterpret_cast<const uint32_t*>(this + 1)[i >> (5 - b)];
3556  //code >>= (i & ((16 >> b) - 1)) << b;
3557  code >>= (i & ((32 >> b) - 1)) << b;
3558  code &= (1 << (1 << b)) - 1;
3559 #endif
3560 #else// use branched version (slow)
3561  float code;
3562  auto *values = reinterpret_cast<const uint8_t*>(this+1);
3563  switch (BaseT::mFlags >> 5) {
3564  case 0u:// 1 bit float
3565  code = float((values[i>>3] >> (i&7) ) & uint8_t(1));
3566  break;
3567  case 1u:// 2 bits float
3568  code = float((values[i>>2] >> ((i&3)<<1)) & uint8_t(3));
3569  break;
3570  case 2u:// 4 bits float
3571  code = float((values[i>>1] >> ((i&1)<<2)) & uint8_t(15));
3572  break;
3573  case 3u:// 8 bits float
3574  code = float(values[i]);
3575  break;
3576  default:// 16 bits float
3577  code = float(reinterpret_cast<const uint16_t*>(values)[i]);
3578  }
3579 #endif
3580  return float(code) * BaseT::mQuantum + BaseT::mMinimum;// code * (max-min)/UNITS + min
3581  }
3582 
3583  /// @brief This class cannot be constructed or deleted
3584  LeafData() = delete;
3585  LeafData(const LeafData&) = delete;
3586  LeafData& operator=(const LeafData&) = delete;
3587  ~LeafData() = delete;
3588 }; // LeafData<FpN>
3589 
3590 // Partial template specialization of LeafData with bool
3591 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3592 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData<bool, CoordT, MaskT, LOG2DIM>
3593 {
3594  static_assert(sizeof(CoordT) == sizeof(Coord), "Mismatching sizeof");
3595  static_assert(sizeof(MaskT<LOG2DIM>) == sizeof(Mask<LOG2DIM>), "Mismatching sizeof");
3596  using ValueType = bool;
3597  using BuildType = bool;
3598  using FloatType = bool;// dummy value type
3599  using ArrayType = MaskT<LOG2DIM>;// type used for the internal mValue array
3600  static constexpr bool FIXED_SIZE = true;
3601 
3602  CoordT mBBoxMin; // 12B.
3603  uint8_t mBBoxDif[3]; // 3B.
3604  uint8_t mFlags; // 1B.
3605  MaskT<LOG2DIM> mValueMask; // LOG2DIM(3): 64B.
3606  MaskT<LOG2DIM> mValues; // LOG2DIM(3): 64B.
3607 
3608  //__hostdev__ const ValueType* values() const { return nullptr; }
3609  __hostdev__ bool getValue(uint32_t i) const { return mValues.isOn(i); }
3610  __hostdev__ bool getMin() const { return false; }// dummy
3611  __hostdev__ bool getMax() const { return false; }// dummy
3612  __hostdev__ bool getAvg() const { return false; }// dummy
3613  __hostdev__ bool getDev() const { return false; }// dummy
3614  __hostdev__ void setValue(uint32_t offset, bool v)
3615  {
3616  mValueMask.setOn(offset);
3617  mValues.set(offset, v);
3618  }
3619 
3620  __hostdev__ void setMin(const bool&) {}// no-op
3621  __hostdev__ void setMax(const bool&) {}// no-op
3622  __hostdev__ void setAvg(const bool&) {}// no-op
3623  __hostdev__ void setDev(const bool&) {}// no-op
3624 
3625  template <typename T>
3626  __hostdev__ void setOrigin(const T& ijk) { mBBoxMin = ijk; }
3627 
3628  /// @brief This class cannot be constructed or deleted
3629  LeafData() = delete;
3630  LeafData(const LeafData&) = delete;
3631  LeafData& operator=(const LeafData&) = delete;
3632  ~LeafData() = delete;
3633 }; // LeafData<bool>
3634 
3635 // Partial template specialization of LeafData with ValueMask
3636 template<typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3637 struct NANOVDB_ALIGN(NANOVDB_DATA_ALIGNMENT) LeafData<ValueMask, CoordT, MaskT, LOG2DIM>
3638 {
3639  static_assert(sizeof(CoordT) == sizeof(Coord), "Mismatching sizeof");
3640  static_assert(sizeof(MaskT<LOG2DIM>) == sizeof(Mask<LOG2DIM>), "Mismatching sizeof");
3641  using ValueType = bool;
3643  using FloatType = bool;// dummy value type
3644  using ArrayType = void;// type used for the internal mValue array - void means missing
3645  static constexpr bool FIXED_SIZE = true;
3646 
3647  CoordT mBBoxMin; // 12B.
3648  uint8_t mBBoxDif[3]; // 3B.
3649  uint8_t mFlags; // 1B.
3650  MaskT<LOG2DIM> mValueMask; // LOG2DIM(3): 64B.
3651 
3652  //__hostdev__ const ValueType* values() const { return nullptr; }
3653  __hostdev__ bool getValue(uint32_t i) const { return mValueMask.isOn(i); }
3654  __hostdev__ bool getMin() const { return false; }// dummy
3655  __hostdev__ bool getMax() const { return false; }// dummy
3656  __hostdev__ bool getAvg() const { return false; }// dummy
3657  __hostdev__ bool getDev() const { return false; }// dummy
3658  __hostdev__ void setValue(uint32_t offset, bool)
3659  {
3660  mValueMask.setOn(offset);
3661  }
3662 
3663  __hostdev__ void setMin(const ValueType&) {}// no-op
3664  __hostdev__ void setMax(const ValueType&) {}// no-op
3665  __hostdev__ void setAvg(const FloatType&) {}// no-op
3666  __hostdev__ void setDev(const FloatType&) {}// no-op
3667 
3668  template <typename T>
3669  __hostdev__ void setOrigin(const T& ijk) { mBBoxMin = ijk; }
3670 
3671  /// @brief This class cannot be constructed or deleted
3672  LeafData() = delete;
3673  LeafData(const LeafData&) = delete;
3674  LeafData& operator=(const LeafData&) = delete;
3675  ~LeafData() = delete;
3676 }; // LeafData<ValueMask>
3677 
3678 /// @brief Leaf nodes of the VDB tree. (defaults to 8x8x8 = 512 voxels)
3679 template<typename BuildT,
3680  typename CoordT = Coord,
3681  template<uint32_t> class MaskT = Mask,
3682  uint32_t Log2Dim = 3>
3683 class LeafNode : private LeafData<BuildT, CoordT, MaskT, Log2Dim>
3684 {
3685 public:
3687  {
3688  __hostdev__ static uint32_t dim() { return 1u; }
3689  }; // Voxel
3692  using ValueType = typename DataType::ValueType;
3693  using FloatType = typename DataType::FloatType;
3694  using BuildType = typename DataType::BuildType;
3695  using CoordType = CoordT;
3696  static constexpr bool FIXED_SIZE = DataType::FIXED_SIZE;
3697  template<uint32_t LOG2>
3698  using MaskType = MaskT<LOG2>;
3699 
3700  static_assert(is_same<ValueType,typename BuildToValueMap<BuildType>::Type>::value, "Mismatching BuildType");
3701  static constexpr uint32_t LOG2DIM = Log2Dim;
3702  static constexpr uint32_t TOTAL = LOG2DIM; // needed by parent nodes
3703  static constexpr uint32_t DIM = 1u << TOTAL; // number of voxels along each axis of this node
3704  static constexpr uint32_t SIZE = 1u << 3 * LOG2DIM; // total number of voxels represented by this node
3705  static constexpr uint32_t MASK = (1u << LOG2DIM) - 1u; // mask for bit operations
3706  static constexpr uint32_t LEVEL = 0; // level 0 = leaf
3707  static constexpr uint64_t NUM_VALUES = uint64_t(1) << (3 * TOTAL); // total voxel count represented by this node
3708 
3709  __hostdev__ DataType* data() { return reinterpret_cast<DataType*>(this); }
3710 
3711  __hostdev__ const DataType* data() const { return reinterpret_cast<const DataType*>(this); }
3712 
3713  /// @brief Return a const reference to the bit mask of active voxels in this leaf node
3714  __hostdev__ const MaskType<LOG2DIM>& valueMask() const { return DataType::mValueMask; }
3715 
3716  /// @brief Return a const reference to the minimum active value encoded in this leaf node
3717  __hostdev__ ValueType minimum() const { return this->getMin(); }
3718 
3719  /// @brief Return a const reference to the maximum active value encoded in this leaf node
3720  __hostdev__ ValueType maximum() const { return this->getMax(); }
3721 
3722  /// @brief Return a const reference to the average of all the active values encoded in this leaf node
3723  __hostdev__ FloatType average() const { return DataType::getAvg(); }
3724 
3725  /// @brief Return the variance of all the active values encoded in this leaf node
3726  __hostdev__ FloatType variance() const { return DataType::getDev()*DataType::getDev(); }
3727 
3728  /// @brief Return a const reference to the standard deviation of all the active values encoded in this leaf node
3729  __hostdev__ FloatType stdDeviation() const { return DataType::getDev(); }
3730 
3731  __hostdev__ uint8_t flags() const { return DataType::mFlags; }
3732 
3733  /// @brief Return the origin in index space of this leaf node
3734  __hostdev__ CoordT origin() const { return DataType::mBBoxMin & ~MASK; }
3735 
3736  __hostdev__ static CoordT OffsetToLocalCoord(uint32_t n)
3737  {
3738  NANOVDB_ASSERT(n < SIZE);
3739  const uint32_t m = n & ((1 << 2 * LOG2DIM) - 1);
3740  return CoordT(n >> 2 * LOG2DIM, m >> LOG2DIM, m & MASK);
3741  }
3742 
3743  /// @brief Converts (in place) a local index coordinate to a global index coordinate
3744  __hostdev__ void localToGlobalCoord(Coord& ijk) const { ijk += this->origin(); }
3745 
3746  __hostdev__ CoordT offsetToGlobalCoord(uint32_t n) const
3747  {
3748  return OffsetToLocalCoord(n) + this->origin();
3749  }
3750 
3751  /// @brief Return the dimension, in index space, of this leaf node (typically 8 as for openvdb leaf nodes!)
3752  __hostdev__ static uint32_t dim() { return 1u << LOG2DIM; }
3753 
3754  /// @brief Return the bounding box in index space of active values in this leaf node
3756  {
3757  BBox<CoordT> bbox(DataType::mBBoxMin, DataType::mBBoxMin);
3758  if ( this->isActive() ) {
3759  bbox.max()[0] += DataType::mBBoxDif[0];
3760  bbox.max()[1] += DataType::mBBoxDif[1];
3761  bbox.max()[2] += DataType::mBBoxDif[2];
3762  } else {// very rare case
3763  bbox = BBox<CoordT>();// invalid
3764  }
3765  return bbox;
3766  }
3767 
3768  /// @brief Return the total number of voxels (e.g. values) encoded in this leaf node
3769  __hostdev__ static uint32_t voxelCount() { return 1u << (3 * LOG2DIM); }
3770 
3771  /// @brief return memory usage in bytes for the class
3772  __hostdev__ static uint64_t memUsage() { return sizeof(LeafNodeType); }
3773 
3774  /// @brief This class cannot be constructed or deleted
3775  LeafNode() = delete;
3776  LeafNode(const LeafNode&) = delete;
3777  LeafNode& operator=(const LeafNode&) = delete;
3778  ~LeafNode() = delete;
3779 
3780  /// @brief Return the voxel value at the given offset.
3781  __hostdev__ ValueType getValue(uint32_t offset) const { return DataType::getValue(offset); }
3782 
3783  /// @brief Return the voxel value at the given coordinate.
3784  __hostdev__ ValueType getValue(const CoordT& ijk) const { return DataType::getValue(CoordToOffset(ijk)); }
3785 
3786  /// @brief Sets the value at the specified location and activate its state.
3787  ///
3788  /// @note This is safe since it does not change the topology of the tree (unlike setValue methods on the other nodes)
3789  __hostdev__ void setValue(const CoordT& ijk, const ValueType& v) { DataType::setValue(CoordToOffset(ijk), v); }
3790 
3791  /// @brief Sets the value at the specified location but leaves its state unchanged.
3792  ///
3793  /// @note This is safe since it does not change the topology of the tree (unlike setValue methods on the other nodes)
3794  __hostdev__ void setValueOnly(uint32_t offset, const ValueType& v) { DataType::setValueOnly(offset, v); }
3795  __hostdev__ void setValueOnly(const CoordT& ijk, const ValueType& v) { DataType::setValueOnly(CoordToOffset(ijk), v); }
3796 
3797  /// @brief Return @c true if the voxel value at the given coordinate is active.
3798  __hostdev__ bool isActive(const CoordT& ijk) const { return DataType::mValueMask.isOn(CoordToOffset(ijk)); }
3799  __hostdev__ bool isActive(uint32_t n) const { return DataType::mValueMask.isOn(n); }
3800 
3801  /// @brief Return @c true if any of the voxel value are active in this leaf node.
3802  __hostdev__ bool isActive() const
3803  {
3804  NANOVDB_ASSERT( bool(DataType::mFlags & uint8_t(2)) != DataType::mValueMask.isOff() );
3805  return DataType::mFlags & uint8_t(2);
3806  }
3807 
3808 
3809  /// @brief Return @c true if the voxel value at the given coordinate is active and updates @c v with the value.
3810  __hostdev__ bool probeValue(const CoordT& ijk, ValueType& v) const
3811  {
3812  const uint32_t n = CoordToOffset(ijk);
3813  v = DataType::getValue(n);
3814  return DataType::mValueMask.isOn(n);
3815  }
3816 
3817  __hostdev__ const LeafNode* probeLeaf(const CoordT&) const { return this; }
3818 
3819  /// @brief Return the linear offset corresponding to the given coordinate
3820  __hostdev__ static uint32_t CoordToOffset(const CoordT& ijk)
3821  {
3822  #if 0
3823  return ((ijk[0] & MASK) << (2 * LOG2DIM)) + ((ijk[1] & MASK) << LOG2DIM) + (ijk[2] & MASK);
3824  #else
3825  return ((ijk[0] & MASK) << (2 * LOG2DIM)) | ((ijk[1] & MASK) << LOG2DIM) | (ijk[2] & MASK);
3826  #endif
3827  }
3828 
3829  /// @brief Updates the local bounding box of active voxels in this node.
3830  ///
3831  /// @warning It assumes that the origin and value mask have already been set.
3832  ///
3833  /// @details This method is based on few (intrinsic) bit operations and hence is relatively fast.
3834  /// However, it should only only be called of either the value mask has changed or if the
3835  /// active bounding box is still undefined. e.g. during constrution of this node.
3836  __hostdev__ void updateBBox();
3837 
3838 private:
3839  static_assert(sizeof(DataType) % NANOVDB_DATA_ALIGNMENT == 0, "sizeof(LeafData) is misaligned");
3840  //static_assert(offsetof(DataType, mValues) % 32 == 0, "LeafData::mValues is misaligned");
3841 
3842  template<typename, int, int, int>
3843  friend class ReadAccessor;
3844 
3845  template<typename>
3846  friend class RootNode;
3847  template<typename, uint32_t>
3848  friend class InternalNode;
3849 
3850  /// @brief Private method to return a voxel value and update a (dummy) ReadAccessor
3851  template<typename AccT>
3852  __hostdev__ ValueType getValueAndCache(const CoordT& ijk, const AccT&) const { return this->getValue(ijk); }
3853 
3854  /// @brief Return the node information.
3855  template<typename AccT>
3856  __hostdev__ typename AccT::NodeInfo getNodeInfoAndCache(const CoordType& /*ijk*/, const AccT& /*acc*/) const {
3857  using NodeInfoT = typename AccT::NodeInfo;
3858  return NodeInfoT{LEVEL, this->dim(), this->minimum(), this->maximum(),
3859  this->average(), this->stdDeviation(), this->bbox()[0], this->bbox()[1]};
3860  }
3861 
3862  template<typename AccT>
3863  __hostdev__ bool isActiveAndCache(const CoordT& ijk, const AccT&) const { return this->isActive(ijk); }
3864 
3865  template<typename AccT>
3866  __hostdev__ bool probeValueAndCache(const CoordT& ijk, ValueType& v, const AccT&) const { return this->probeValue(ijk, v); }
3867 
3868  template<typename AccT>
3869  __hostdev__ const LeafNode* probeLeafAndCache(const CoordT&, const AccT&) const { return this; }
3870 
3871  template<typename RayT, typename AccT>
3872  __hostdev__ uint32_t getDimAndCache(const CoordT&, const RayT& /*ray*/, const AccT&) const
3873  {
3874  if (DataType::mFlags & uint8_t(1))
3875  return this->dim(); // skip this node if first bit is set
3876  //if (!ray.intersects( this->bbox() )) return 1 << LOG2DIM;
3877  return ChildNodeType::dim();
3878  }
3879 
3880 }; // LeafNode class
3881 
3882 template<typename ValueT, typename CoordT, template<uint32_t> class MaskT, uint32_t LOG2DIM>
3884 {
3885  static_assert(LOG2DIM == 3, "LeafNode::updateBBox: only supports LOGDIM = 3!");
3886  if (!this->isActive()) return;
3887  auto update = [&](uint32_t min, uint32_t max, int axis) {
3888  NANOVDB_ASSERT(min <= max && max < 8);
3889  DataType::mBBoxMin[axis] = (DataType::mBBoxMin[axis] & ~MASK) + int(min);
3890  DataType::mBBoxDif[axis] = uint8_t(max - min);
3891  };
3892  uint64_t word64 = DataType::mValueMask.template getWord<uint64_t>(0);
3893  uint32_t Xmin = word64 ? 0u : 8u;
3894  uint32_t Xmax = Xmin;
3895  for (int i = 1; i < 8; ++i) { // last loop over 8 64 words
3896  if (uint64_t w = DataType::mValueMask.template getWord<uint64_t>(i)) { // skip if word has no set bits
3897  word64 |= w; // union 8 x 64 bits words into one 64 bit word
3898  if (Xmin == 8) {
3899  Xmin = i; // only set once
3900  }
3901  Xmax = i;
3902  }
3903  }
3904  NANOVDB_ASSERT(word64);
3905  update(Xmin, Xmax, 0);
3906  update(FindLowestOn(word64) >> 3, FindHighestOn(word64) >> 3, 1);
3907  const uint32_t *p = reinterpret_cast<const uint32_t*>(&word64), word32 = p[0] | p[1];
3908  const uint16_t *q = reinterpret_cast<const uint16_t*>(&word32), word16 = q[0] | q[1];
3909  const uint8_t *b = reinterpret_cast<const uint8_t* >(&word16), byte = b[0] | b[1];
3910  NANOVDB_ASSERT(byte);
3911  update(FindLowestOn(static_cast<uint32_t>(byte)), FindHighestOn(static_cast<uint32_t>(byte)), 2);
3912 } // LeafNode::updateBBox
3913 
3914 // --------------------------> Template specializations and traits <------------------------------------
3915 
3916 /// @brief Template specializations to the default configuration used in OpenVDB:
3917 /// Root -> 32^3 -> 16^3 -> 8^3
3918 template<typename BuildT>
3920 template<typename BuildT>
3922 template<typename BuildT>
3924 template<typename BuildT>
3926 template<typename BuildT>
3928 template<typename BuildT>
3930 
3931 /// @brief Trait to map from LEVEL to node type
3932 template<typename BuildT, int LEVEL>
3933 struct NanoNode;
3934 
3935 // Partial template specialization of above Node struct
3936 template<typename BuildT>
3937 struct NanoNode<BuildT, 0>
3938 {
3941 };
3942 template<typename BuildT>
3943 struct NanoNode<BuildT, 1>
3944 {
3947 };
3948 template<typename BuildT>
3949 struct NanoNode<BuildT, 2>
3950 {
3953 };
3954 template<typename BuildT>
3955 struct NanoNode<BuildT, 3>
3956 {
3959 };
3960 
3973 
3986 
3987 // --------------------------> ReadAccessor <------------------------------------
3988 
3989 /// @brief A read-only value accessor with three levels of node caching. This allows for
3990 /// inverse tree traversal during lookup, which is on average significantly faster
3991 /// than calling the equivalent method on the tree (i.e. top-down traversal).
3992 ///
3993 /// @note By virtue of the fact that a value accessor accelerates random access operations
3994 /// by re-using cached access patterns, this access should be reused for multiple access
3995 /// operations. In other words, never create an instance of this accessor for a single
3996 /// acccess only. In general avoid single access operations with this accessor, and
3997 /// if that is not possible call the corresponding method on the tree instead.
3998 ///
3999 /// @warning Since this ReadAccessor internally caches raw pointers to the nodes of the tree
4000 /// structure, it is not safe to copy between host and device, or even to share among
4001 /// multiple threads on the same host or device. However, it is light-weight so simple
4002 /// instantiate one per thread (on the host and/or device).
4003 ///
4004 /// @details Used to accelerated random access into a VDB tree. Provides on average
4005 /// O(1) random access operations by means of inverse tree traversal,
4006 /// which amortizes the non-const time complexity of the root node.
4007 
4008 template <typename BuildT>
4009 class ReadAccessor<BuildT, -1, -1, -1>
4010 {
4011  using RootT = NanoRoot<BuildT>; // root node
4012  using LeafT = NanoLeaf<BuildT>; // Leaf node
4013  using FloatType = typename RootT::FloatType;
4014  using CoordValueType = typename RootT::CoordType::ValueType;
4015 
4016  mutable const RootT* mRoot; // 8 bytes (mutable to allow for access methods to be const)
4017 public:
4018  using ValueType = typename RootT::ValueType;
4019  using CoordType = typename RootT::CoordType;
4020 
4021  static const int CacheLevels = 0;
4022 
4023  struct NodeInfo {
4024  uint32_t mLevel; // 4B
4025  uint32_t mDim; // 4B
4026  ValueType mMinimum; // typically 4B
4027  ValueType mMaximum; // typically 4B
4028  FloatType mAverage; // typically 4B
4029  FloatType mStdDevi; // typically 4B
4032  };
4033 
4034  /// @brief Constructor from a root node
4035  __hostdev__ ReadAccessor(const RootT& root) : mRoot{&root} {}
4036 
4037  __hostdev__ const RootT& root() const { return *mRoot; }
4038 
4039  /// @brief Defaults constructors
4040  ReadAccessor(const ReadAccessor&) = default;
4041  ~ReadAccessor() = default;
4042  ReadAccessor& operator=(const ReadAccessor&) = default;
4043 
4045  {
4046  return mRoot->getValueAndCache(ijk, *this);
4047  }
4048 
4049  __hostdev__ NodeInfo getNodeInfo(const CoordType& ijk) const
4050  {
4051  return mRoot->getNodeInfoAndCache(ijk, *this);
4052  }
4053 
4054  __hostdev__ bool isActive(const CoordType& ijk) const
4055  {
4056  return mRoot->isActiveAndCache(ijk, *this);
4057  }
4058 
4059  __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const
4060  {
4061  return mRoot->probeValueAndCache(ijk, v, *this);
4062  }
4063 
4064  __hostdev__ const LeafT* probeLeaf(const CoordType& ijk) const
4065  {
4066  return mRoot->probeLeafAndCache(ijk, *this);
4067  }
4068 
4069  template<typename RayT>
4070  __hostdev__ uint32_t getDim(const CoordType& ijk, const RayT& ray) const
4071  {
4072  return mRoot->getDimAndCache(ijk, ray, *this);
4073  }
4074 
4075 private:
4076  /// @brief Allow nodes to insert themselves into the cache.
4077  template<typename>
4078  friend class RootNode;
4079  template<typename, uint32_t>
4080  friend class InternalNode;
4081  template<typename, typename, template<uint32_t> class, uint32_t>
4082  friend class LeafNode;
4083 
4084  /// @brief No-op
4085  template<typename NodeT>
4086  __hostdev__ void insert(const CoordType&, const NodeT*) const {}
4087 }; // ReadAccessor<ValueT, -1, -1, -1> class
4088 
4089 /// @brief Node caching at a single tree level
4090 template <typename BuildT, int LEVEL0>
4091 class ReadAccessor<BuildT, LEVEL0, -1, -1>//e.g. 0, 1, 2
4092 {
4093  static_assert(LEVEL0 >= 0 && LEVEL0 <= 2, "LEVEL0 should be 0, 1, or 2");
4094 
4095  using TreeT = NanoTree<BuildT>;
4096  using RootT = NanoRoot<BuildT>; // root node
4097  using LeafT = NanoLeaf<BuildT>; // Leaf node
4098  using NodeT = typename NodeTrait<TreeT, LEVEL0>::type;
4099  using CoordT = typename RootT::CoordType;
4100  using ValueT = typename RootT::ValueType;
4101 
4102  using FloatType = typename RootT::FloatType;
4103  using CoordValueType = typename RootT::CoordT::ValueType;
4104 
4105  // All member data are mutable to allow for access methods to be const
4106  mutable CoordT mKey; // 3*4 = 12 bytes
4107  mutable const RootT* mRoot; // 8 bytes
4108  mutable const NodeT* mNode; // 8 bytes
4109 
4110 public:
4111  using ValueType = ValueT;
4112  using CoordType = CoordT;
4113 
4114  static const int CacheLevels = 1;
4115 
4116  using NodeInfo = typename ReadAccessor<ValueT, -1, -1, -1>::NodeInfo;
4117 
4118  /// @brief Constructor from a root node
4120  : mKey(CoordType::max())
4121  , mRoot(&root)
4122  , mNode(nullptr)
4123  {
4124  }
4125 
4126  __hostdev__ const RootT& root() const { return *mRoot; }
4127 
4128  /// @brief Defaults constructors
4129  ReadAccessor(const ReadAccessor&) = default;
4130  ~ReadAccessor() = default;
4131  ReadAccessor& operator=(const ReadAccessor&) = default;
4132 
4133  __hostdev__ bool isCached(const CoordType& ijk) const
4134  {
4135  return (ijk[0] & int32_t(~NodeT::MASK)) == mKey[0] &&
4136  (ijk[1] & int32_t(~NodeT::MASK)) == mKey[1] &&
4137  (ijk[2] & int32_t(~NodeT::MASK)) == mKey[2];
4138  }
4139 
4141  {
4142  if (this->isCached(ijk)) {
4143  return mNode->getValueAndCache(ijk, *this);
4144  }
4145  return mRoot->getValueAndCache(ijk, *this);
4146  }
4147 
4149  {
4150  if (this->isCached(ijk)) {
4151  return mNode->getNodeInfoAndCache(ijk, *this);
4152  }
4153  return mRoot->getNodeInfoAndCache(ijk, *this);
4154  }
4155 
4156  __hostdev__ bool isActive(const CoordType& ijk) const
4157  {
4158  if (this->isCached(ijk)) {
4159  return mNode->isActiveAndCache(ijk, *this);
4160  }
4161  return mRoot->isActiveAndCache(ijk, *this);
4162  }
4163 
4164  __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const
4165  {
4166  if (this->isCached(ijk)) {
4167  return mNode->probeValueAndCache(ijk, v, *this);
4168  }
4169  return mRoot->probeValueAndCache(ijk, v, *this);
4170  }
4171 
4172  __hostdev__ const LeafT* probeLeaf(const CoordType& ijk) const
4173  {
4174  if (this->isCached(ijk)) {
4175  return mNode->probeLeafAndCache(ijk, *this);
4176  }
4177  return mRoot->probeLeafAndCache(ijk, *this);
4178  }
4179 
4180  template<typename RayT>
4181  __hostdev__ uint32_t getDim(const CoordType& ijk, const RayT& ray) const
4182  {
4183  if (this->isCached(ijk)) {
4184  return mNode->getDimAndCache(ijk, ray, *this);
4185  }
4186  return mRoot->getDimAndCache(ijk, ray, *this);
4187  }
4188 
4189 private:
4190  /// @brief Allow nodes to insert themselves into the cache.
4191  template<typename>
4192  friend class RootNode;
4193  template<typename, uint32_t>
4194  friend class InternalNode;
4195  template<typename, typename, template<uint32_t> class, uint32_t>
4196  friend class LeafNode;
4197 
4198  /// @brief Inserts a leaf node and key pair into this ReadAccessor
4199  __hostdev__ void insert(const CoordType& ijk, const NodeT* node) const
4200  {
4201  mKey = ijk & ~NodeT::MASK;
4202  mNode = node;
4203  }
4204 
4205  // no-op
4206  template<typename OtherNodeT>
4207  __hostdev__ void insert(const CoordType&, const OtherNodeT*) const {}
4208 
4209 }; // ReadAccessor<ValueT, LEVEL0>
4210 
4211 template <typename BuildT, int LEVEL0, int LEVEL1>
4212 class ReadAccessor<BuildT, LEVEL0, LEVEL1, -1>//e.g. (0,1), (1,2), (0,2)
4213 {
4214  static_assert(LEVEL0 >= 0 && LEVEL0 <= 2, "LEVEL0 must be 0, 1, 2");
4215  static_assert(LEVEL1 >= 0 && LEVEL1 <= 2, "LEVEL1 must be 0, 1, 2");
4216  static_assert(LEVEL0 < LEVEL1, "Level 0 must be lower than level 1");
4217  using TreeT = NanoTree<BuildT>;
4218  using RootT = NanoRoot<BuildT>;
4219  using LeafT = NanoLeaf<BuildT>;
4220  using Node1T = typename NodeTrait<TreeT, LEVEL0>::type;
4221  using Node2T = typename NodeTrait<TreeT, LEVEL1>::type;
4222  using CoordT = typename RootT::CoordType;
4223  using ValueT = typename RootT::ValueType;
4224  using FloatType = typename RootT::FloatType;
4225  using CoordValueType = typename RootT::CoordT::ValueType;
4226 
4227  // All member data are mutable to allow for access methods to be const
4228 #ifdef USE_SINGLE_ACCESSOR_KEY // 44 bytes total
4229  mutable CoordT mKey; // 3*4 = 12 bytes
4230 #else // 68 bytes total
4231  mutable CoordT mKeys[2]; // 2*3*4 = 24 bytes
4232 #endif
4233  mutable const RootT* mRoot;
4234  mutable const Node1T* mNode1;
4235  mutable const Node2T* mNode2;
4236 
4237 public:
4238  using ValueType = ValueT;
4239  using CoordType = CoordT;
4240 
4241  static const int CacheLevels = 2;
4242 
4243  using NodeInfo = typename ReadAccessor<ValueT,-1,-1,-1>::NodeInfo;
4244 
4245  /// @brief Constructor from a root node
4247 #ifdef USE_SINGLE_ACCESSOR_KEY
4248  : mKey(CoordType::max())
4249 #else
4250  : mKeys{CoordType::max(), CoordType::max()}
4251 #endif
4252  , mRoot(&root)
4253  , mNode1(nullptr)
4254  , mNode2(nullptr)
4255  {
4256  }
4257 
4258  __hostdev__ const RootT& root() const { return *mRoot; }
4259 
4260  /// @brief Defaults constructors
4261  ReadAccessor(const ReadAccessor&) = default;
4262  ~ReadAccessor() = default;
4263  ReadAccessor& operator=(const ReadAccessor&) = default;
4264 
4265 #ifdef USE_SINGLE_ACCESSOR_KEY
4266  __hostdev__ bool isCached1(CoordValueType dirty) const
4267  {
4268  if (!mNode1)
4269  return false;
4270  if (dirty & int32_t(~Node1T::MASK)) {
4271  mNode1 = nullptr;
4272  return false;
4273  }
4274  return true;
4275  }
4276  __hostdev__ bool isCached2(CoordValueType dirty) const
4277  {
4278  if (!mNode2)
4279  return false;
4280  if (dirty & int32_t(~Node2T::MASK)) {
4281  mNode2 = nullptr;
4282  return false;
4283  }
4284  return true;
4285  }
4286  __hostdev__ CoordValueType computeDirty(const CoordType& ijk) const
4287  {
4288  return (ijk[0] ^ mKey[0]) | (ijk[1] ^ mKey[1]) | (ijk[2] ^ mKey[2]);
4289  }
4290 #else
4291  __hostdev__ bool isCached1(const CoordType& ijk) const
4292  {
4293  return (ijk[0] & int32_t(~Node1T::MASK)) == mKeys[0][0] &&
4294  (ijk[1] & int32_t(~Node1T::MASK)) == mKeys[0][1] &&
4295  (ijk[2] & int32_t(~Node1T::MASK)) == mKeys[0][2];
4296  }
4297  __hostdev__ bool isCached2(const CoordType& ijk) const
4298  {
4299  return (ijk[0] & int32_t(~Node2T::MASK)) == mKeys[1][0] &&
4300  (ijk[1] & int32_t(~Node2T::MASK)) == mKeys[1][1] &&
4301  (ijk[2] & int32_t(~Node2T::MASK)) == mKeys[1][2];
4302  }
4303 #endif
4304 
4306  {
4307 #ifdef USE_SINGLE_ACCESSOR_KEY
4308  const CoordValueType dirty = this->computeDirty(ijk);
4309 #else
4310  auto&& dirty = ijk;
4311 #endif
4312  if (this->isCached1(dirty)) {
4313  return mNode1->getValueAndCache(ijk, *this);
4314  } else if (this->isCached2(dirty)) {
4315  return mNode2->getValueAndCache(ijk, *this);
4316  }
4317  return mRoot->getValueAndCache(ijk, *this);
4318  }
4319 
4321  {
4322 #ifdef USE_SINGLE_ACCESSOR_KEY
4323  const CoordValueType dirty = this->computeDirty(ijk);
4324 #else
4325  auto&& dirty = ijk;
4326 #endif
4327  if (this->isCached1(dirty)) {
4328  return mNode1->getNodeInfoAndCache(ijk, *this);
4329  } else if (this->isCached2(dirty)) {
4330  return mNode2->getNodeInfoAndCache(ijk, *this);
4331  }
4332  return mRoot->getNodeInfoAndCache(ijk, *this);
4333  }
4334 
4335  __hostdev__ bool isActive(const CoordType& ijk) const
4336  {
4337 #ifdef USE_SINGLE_ACCESSOR_KEY
4338  const CoordValueType dirty = this->computeDirty(ijk);
4339 #else
4340  auto&& dirty = ijk;
4341 #endif
4342  if (this->isCached1(dirty)) {
4343  return mNode1->isActiveAndCache(ijk, *this);
4344  } else if (this->isCached2(dirty)) {
4345  return mNode2->isActiveAndCache(ijk, *this);
4346  }
4347  return mRoot->isActiveAndCache(ijk, *this);
4348  }
4349 
4350  __hostdev__ bool probeValue(const CoordType& ijk, ValueType& v) const
4351  {
4352 #ifdef USE_SINGLE_ACCESSOR_KEY
4353  const CoordValueType dirty = this->computeDirty(ijk);
4354 #else
4355  auto&& dirty = ijk;
4356 #endif
4357  if (this->isCached1(dirty)) {
4358  return mNode1->probeValueAndCache(ijk, v, *this);
4359  } else if (this->isCached2(dirty)) {
4360  return mNode2->probeValueAndCache(ijk, v, *this);
4361  }
4362  return mRoot->probeValueAndCache(ijk, v, *this);
4363  }
4364 
4365  __hostdev__ const LeafT* probeLeaf(const CoordType& ijk) const
4366  {
4367 #ifdef USE_SINGLE_ACCESSOR_KEY
4368  const CoordValueType dirty = this->computeDirty(ijk);
4369 #else
4370  auto&& dirty = ijk;
4371 #endif
4372  if (this->isCached1(dirty)) {
4373  return mNode1->probeLeafAndCache(ijk, *this);
4374  } else if (this->isCached2(dirty)) {
4375  return mNode2->probeLeafAndCache(ijk, *this);
4376  }
4377  return mRoot->probeLeafAndCache(ijk, *this);
4378  }
4379 
4380  template<typename RayT>
4381  __hostdev__ uint32_t getDim(const CoordType& ijk, const RayT& ray) const
4382  {
4383 #ifdef USE_SINGLE_ACCESSOR_KEY
4384  const CoordValueType dirty = this->computeDirty(ijk);
4385 #else
4386  auto&& dirty = ijk;
4387 #endif
4388  if (this->isCached1(dirty)) {
4389  return mNode1->getDimAndCache(ijk, ray, *this);
4390  } else if (this->isCached2(dirty)) {
4391  return mNode2->getDimAndCache(ijk, ray, *this);
4392  }
4393  return mRoot->getDimAndCache(ijk, ray, *this);
4394  }
4395 
4396 private:
4397  /// @brief Allow nodes to insert themselves into the cache.
4398  template<typename>
4399  friend class RootNode;
4400  template<typename, uint32_t>
4401  friend class InternalNode;
4402  template<typename, typename, template<uint32_t> class, uint32_t>
4403  friend class LeafNode;
4404 
4405  /// @brief Inserts a leaf node and key pair into this ReadAccessor
4406  __hostdev__ void insert(const CoordType& ijk, const Node1T* node) const
4407  {
4408 #ifdef USE_SINGLE_ACCESSOR_KEY
4409  mKey = ijk;
4410 #else
4411  mKeys[0] = ijk & ~Node1T::MASK;
4412 #endif
4413  mNode1 = node;
4414  }
4415  __hostdev__ void insert(const CoordType& ijk, const Node2T* node) const
4416  {
4417 #ifdef USE_SINGLE_ACCESSOR_KEY
4418  mKey = ijk;
4419 #else
4420  mKeys[1] = ijk & ~Node2T::MASK;
4421 #endif
4422  mNode2 = node;
4423  }
4424  template <typename OtherNodeT>
4425  __hostdev__ void insert(const CoordType&, const OtherNodeT*) const {}
4426 }; // ReadAccessor<BuildT, LEVEL0, LEVEL1>
4427 
4428 
4429 /// @brief Node caching at all (three) tree levels
4430 template <typename BuildT>
4431 class ReadAccessor<BuildT, 0, 1, 2>
4432 {
4433  using TreeT = NanoTree<BuildT>;
4434  using RootT = NanoRoot<BuildT>; // root node
4435  using NodeT2 = NanoUpper<BuildT>; // upper internal node
4436  using NodeT1 = NanoLower<BuildT>; // lower internal node
4437  using LeafT = NanoLeaf< BuildT>; // Leaf node
4438  using CoordT = typename RootT::CoordType;
4439  using ValueT = typename RootT::ValueType;
4440 
4441  using FloatType = typename RootT::FloatType;
4442  using CoordValueType = typename RootT::CoordT::ValueType;
4443 
4444  // All member data are mutable to allow for access methods to be const
4445 #ifdef USE_SINGLE_ACCESSOR_KEY // 44 bytes total
4446  mutable CoordT mKey; // 3*4 = 12 bytes
4447 #else // 68 bytes total
4448  mutable CoordT mKeys[3]; // 3*3*4 = 36 bytes
4449 #endif
4450  mutable const RootT* mRoot;
4451  mutable const void* mNode[3]; // 4*8 = 32 bytes
4452 
4453 public:
4454  using ValueType = ValueT;
4455  using CoordType = CoordT;
4456 
4457  static const int CacheLevels = 3;
4458 
4459  using NodeInfo = typename ReadAccessor<ValueT, -1, -1, -1>::NodeInfo;
4460 
4461  /// @brief Constructor from a root node
4463 #ifdef USE_SINGLE_ACCESSOR_KEY
4464  : mKey(CoordType::max())
4465 #else
4467 #endif
4468  , mRoot(&root)
4469  , mNode{nullptr, nullptr, nullptr}
4470  {
4471  }
4472 
4473  __hostdev__ const RootT& root() const { return *mRoot; }
4474 
4475  /// @brief Defaults constructors
4476  ReadAccessor(const ReadAccessor&) = default;
4477  ~ReadAccessor() = default;
4478  ReadAccessor& operator=(const ReadAccessor&) = default;
4479 
4480  /// @brief Return a const point to the cached node of the specified type
4481  ///
4482  /// @warning The return value could be NULL.
4483  template<typename NodeT>
4484  __hostdev__ const NodeT* getNode() const
4485  {
4486  using T = typename No