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