HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SampleFromVoxels.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 ///
6 /// @file SampleFromVoxels.h
7 ///
8 /// @brief NearestNeighborSampler, TrilinearSampler, TriquadraticSampler and TricubicSampler
9 ///
10 /// @note These interpolators employ internal caching for better performance when used repeatedly
11 /// in the same voxel location, so try to reuse an instance of these classes more than once.
12 ///
13 /// @warning While all the interpolators defined below work with both scalars and vectors
14 /// values (e.g. float and Vec3<float>) TrilinarSampler::zeroCrossing and
15 /// Trilinear::gradient will only compile with floating point value types.
16 ///
17 /// @author Ken Museth
18 ///
19 ///////////////////////////////////////////////////////////////////////////
20 
21 #ifndef NANOVDB_SAMPLE_FROM_VOXELS_H_HAS_BEEN_INCLUDED
22 #define NANOVDB_SAMPLE_FROM_VOXELS_H_HAS_BEEN_INCLUDED
23 
24 // Only define __hostdev__ when compiling as NVIDIA CUDA
25 #if defined(__CUDACC__) || defined(__HIP__)
26 #define __hostdev__ __host__ __device__
27 #else
28 #include <cmath> // for floor
29 #define __hostdev__
30 #endif
31 
32 namespace nanovdb {
33 
34 // Forward declaration of sampler with specific polynomial orders
35 template<typename TreeT, int Order, bool UseCache = true>
36 class SampleFromVoxels;
37 
38 /// @brief Factory free-function for a sampler of specific polynomial orders
39 ///
40 /// @details This allows for the compact syntax:
41 /// @code
42 /// auto acc = grid.getAccessor();
43 /// auto smp = nanovdb::createSampler<1>( acc );
44 /// @endcode
45 template<int Order, typename TreeOrAccT, bool UseCache = true>
46 __hostdev__ SampleFromVoxels<TreeOrAccT, Order, UseCache> createSampler(const TreeOrAccT& acc)
47 {
48  return SampleFromVoxels<TreeOrAccT, Order, UseCache>(acc);
49 }
50 
51 /// @brief Utility function that returns the Coord of the round-down of @a xyz
52 /// and redefined @xyz as the fractional part, ie xyz-in = return-value + xyz-out
53 template<typename CoordT, typename RealT, template<typename> class Vec3T>
54 __hostdev__ inline CoordT Floor(Vec3T<RealT>& xyz);
55 
56 /// @brief Template specialization of Floor for Vec3<float>
57 template<typename CoordT, template<typename> class Vec3T>
58 __hostdev__ inline CoordT Floor(Vec3T<float>& xyz)
59 {
60  const float ijk[3] = {floorf(xyz[0]), floorf(xyz[1]), floorf(xyz[2])};
61  xyz[0] -= ijk[0];
62  xyz[1] -= ijk[1];
63  xyz[2] -= ijk[2];
64  return CoordT(int32_t(ijk[0]), int32_t(ijk[1]), int32_t(ijk[2]));
65 }
66 
67 /// @brief Template specialization of Floor for Vec3<float>
68 template<typename CoordT, template<typename> class Vec3T>
69 __hostdev__ inline CoordT Floor(Vec3T<double>& xyz)
70 {
71  const double ijk[3] = {floor(xyz[0]), floor(xyz[1]), floor(xyz[2])};
72  xyz[0] -= ijk[0];
73  xyz[1] -= ijk[1];
74  xyz[2] -= ijk[2];
75  return CoordT(int32_t(ijk[0]), int32_t(ijk[1]), int32_t(ijk[2]));
76 }
77 
78 // ------------------------------> NearestNeighborSampler <--------------------------------------
79 
80 /// @brief Nearest neighbor, i.e. zero order, interpolator with caching
81 template<typename TreeOrAccT>
82 class SampleFromVoxels<TreeOrAccT, 0, true>
83 {
84 public:
85  using ValueT = typename TreeOrAccT::ValueType;
86  using CoordT = typename TreeOrAccT::CoordType;
87 
88  static const int ORDER = 0;
89  /// @brief Construction from a Tree or ReadAccessor
90  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc)
91  : mAcc(acc)
92  , mPos(CoordT::max())
93  {
94  }
95 
96  __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
97 
98  /// @note xyz is in index space space
99  template<typename Vec3T>
100  inline __hostdev__ ValueT operator()(const Vec3T& xyz) const;
101 
102  inline __hostdev__ ValueT operator()(const CoordT& ijk) const;
103 
104 private:
105  const TreeOrAccT& mAcc;
106  mutable CoordT mPos;
107  mutable ValueT mVal; // private cache
108 }; // SampleFromVoxels<TreeOrAccT, 0, true>
109 
110 /// @brief Nearest neighbor, i.e. zero order, interpolator without caching
111 template<typename TreeOrAccT>
112 class SampleFromVoxels<TreeOrAccT, 0, false>
113 {
114 public:
115  using ValueT = typename TreeOrAccT::ValueType;
116  using CoordT = typename TreeOrAccT::CoordType;
117  static const int ORDER = 0;
118 
119  /// @brief Construction from a Tree or ReadAccessor
120  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc)
121  : mAcc(acc)
122  {
123  }
124 
125  __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
126 
127  /// @note xyz is in index space space
128  template<typename Vec3T>
129  inline __hostdev__ ValueT operator()(const Vec3T& xyz) const;
130 
131  inline __hostdev__ ValueT operator()(const CoordT& ijk) const { return mAcc.getValue(ijk);}
132 
133 private:
134  const TreeOrAccT& mAcc;
135 }; // SampleFromVoxels<TreeOrAccT, 0, false>
136 
137 template<typename TreeOrAccT>
138 template<typename Vec3T>
139 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 0, true>::operator()(const Vec3T& xyz) const
140 {
141  const CoordT ijk = Round<CoordT>(xyz);
142  if (ijk != mPos) {
143  mPos = ijk;
144  mVal = mAcc.getValue(mPos);
145  }
146  return mVal;
147 }
148 
149 template<typename TreeOrAccT>
150 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 0, true>::operator()(const CoordT& ijk) const
151 {
152  if (ijk != mPos) {
153  mPos = ijk;
154  mVal = mAcc.getValue(mPos);
155  }
156  return mVal;
157 }
158 
159 template<typename TreeOrAccT>
160 template<typename Vec3T>
161 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 0, false>::operator()(const Vec3T& xyz) const
162 {
163  return mAcc.getValue(Round<CoordT>(xyz));
164 }
165 
166 // ------------------------------> TrilinearSampler <--------------------------------------
167 
168 /// @brief Tri-linear sampler, i.e. first order, interpolator
169 template<typename TreeOrAccT>
170 class TrilinearSampler
171 {
172 protected:
173  const TreeOrAccT& mAcc;
174 
175 public:
176  using ValueT = typename TreeOrAccT::ValueType;
177  using CoordT = typename TreeOrAccT::CoordType;
178  static const int ORDER = 1;
179 
180  /// @brief Protected constructor from a Tree or ReadAccessor
181  __hostdev__ TrilinearSampler(const TreeOrAccT& acc) : mAcc(acc) {}
182 
183  __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
184 
185  /// @brief Extract the stencil of 8 values
186  inline __hostdev__ void stencil(CoordT& ijk, ValueT (&v)[2][2][2]) const;
187 
188  template<typename RealT, template<typename...> class Vec3T>
189  static inline __hostdev__ ValueT sample(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2]);
190 
191  template<typename RealT, template<typename...> class Vec3T>
192  static inline __hostdev__ Vec3T<ValueT> gradient(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2]);
193 
194  static inline __hostdev__ bool zeroCrossing(const ValueT (&v)[2][2][2]);
195 }; // TrilinearSamplerBase
196 
197 template<typename TreeOrAccT>
198 __hostdev__ void TrilinearSampler<TreeOrAccT>::stencil(CoordT& ijk, ValueT (&v)[2][2][2]) const
199 {
200  v[0][0][0] = mAcc.getValue(ijk); // i, j, k
201 
202  ijk[2] += 1;
203  v[0][0][1] = mAcc.getValue(ijk); // i, j, k + 1
204 
205  ijk[1] += 1;
206  v[0][1][1] = mAcc.getValue(ijk); // i, j+1, k + 1
207 
208  ijk[2] -= 1;
209  v[0][1][0] = mAcc.getValue(ijk); // i, j+1, k
210 
211  ijk[0] += 1;
212  ijk[1] -= 1;
213  v[1][0][0] = mAcc.getValue(ijk); // i+1, j, k
214 
215  ijk[2] += 1;
216  v[1][0][1] = mAcc.getValue(ijk); // i+1, j, k + 1
217 
218  ijk[1] += 1;
219  v[1][1][1] = mAcc.getValue(ijk); // i+1, j+1, k + 1
220 
221  ijk[2] -= 1;
222  v[1][1][0] = mAcc.getValue(ijk); // i+1, j+1, k
223 }
224 
225 template<typename TreeOrAccT>
226 template<typename RealT, template<typename...> class Vec3T>
227 __hostdev__ typename TreeOrAccT::ValueType TrilinearSampler<TreeOrAccT>::sample(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2])
228 {
229 #if 0
230  auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b-a, a); };// = w*(b-a) + a
231  //auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b, fma(-w, a, a));};// = (1-w)*a + w*b
232 #else
233  auto lerp = [](ValueT a, ValueT b, RealT w) { return a + ValueT(w) * (b - a); };
234 #endif
235  return lerp(lerp(lerp(v[0][0][0], v[0][0][1], uvw[2]), lerp(v[0][1][0], v[0][1][1], uvw[2]), uvw[1]),
236  lerp(lerp(v[1][0][0], v[1][0][1], uvw[2]), lerp(v[1][1][0], v[1][1][1], uvw[2]), uvw[1]),
237  uvw[0]);
238 }
239 
240 template<typename TreeOrAccT>
241 template<typename RealT, template<typename...> class Vec3T>
242 __hostdev__ Vec3T<typename TreeOrAccT::ValueType> TrilinearSampler<TreeOrAccT>::gradient(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2])
243 {
244  static_assert(is_floating_point<ValueT>::value, "TrilinearSampler::gradient requires a floating-point type");
245 #if 0
246  auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b-a, a); };// = w*(b-a) + a
247  //auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b, fma(-w, a, a));};// = (1-w)*a + w*b
248 #else
249  auto lerp = [](ValueT a, ValueT b, RealT w) { return a + ValueT(w) * (b - a); };
250 #endif
251 
252  ValueT D[4] = {v[0][0][1] - v[0][0][0], v[0][1][1] - v[0][1][0], v[1][0][1] - v[1][0][0], v[1][1][1] - v[1][1][0]};
253 
254  // Z component
255  Vec3T<ValueT> grad(0, 0, lerp(lerp(D[0], D[1], uvw[1]), lerp(D[2], D[3], uvw[1]), uvw[0]));
256 
257  const ValueT w = ValueT(uvw[2]);
258  D[0] = v[0][0][0] + D[0] * w;
259  D[1] = v[0][1][0] + D[1] * w;
260  D[2] = v[1][0][0] + D[2] * w;
261  D[3] = v[1][1][0] + D[3] * w;
262 
263  // X component
264  grad[0] = lerp(D[2], D[3], uvw[1]) - lerp(D[0], D[1], uvw[1]);
265 
266  // Y component
267  grad[1] = lerp(D[1] - D[0], D[3] - D[2], uvw[0]);
268 
269  return grad;
270 }
271 
272 template<typename TreeOrAccT>
273 __hostdev__ bool TrilinearSampler<TreeOrAccT>::zeroCrossing(const ValueT (&v)[2][2][2])
274 {
275  static_assert(is_floating_point<ValueT>::value, "TrilinearSampler::zeroCrossing requires a floating-point type");
276  const bool less = v[0][0][0] < ValueT(0);
277  return (less ^ (v[0][0][1] < ValueT(0))) ||
278  (less ^ (v[0][1][1] < ValueT(0))) ||
279  (less ^ (v[0][1][0] < ValueT(0))) ||
280  (less ^ (v[1][0][0] < ValueT(0))) ||
281  (less ^ (v[1][0][1] < ValueT(0))) ||
282  (less ^ (v[1][1][1] < ValueT(0))) ||
283  (less ^ (v[1][1][0] < ValueT(0)));
284 }
285 
286 /// @brief Template specialization that does not use caching of stencil points
287 template<typename TreeOrAccT>
288 class SampleFromVoxels<TreeOrAccT, 1, false> : public TrilinearSampler<TreeOrAccT>
289 {
290  using BaseT = TrilinearSampler<TreeOrAccT>;
291  using ValueT = typename TreeOrAccT::ValueType;
292  using CoordT = typename TreeOrAccT::CoordType;
293 
294 public:
295 
296  /// @brief Construction from a Tree or ReadAccessor
297  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc) {}
298 
299  /// @note xyz is in index space space
300  template<typename RealT, template<typename...> class Vec3T>
301  inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
302 
303  /// @note ijk is in index space space
304  __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
305 
306  /// @brief Return the gradient in index space.
307  ///
308  /// @warning Will only compile with floating point value types
309  template<typename RealT, template<typename...> class Vec3T>
310  inline __hostdev__ Vec3T<ValueT> gradient(Vec3T<RealT> xyz) const;
311 
312  /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
313  ///
314  /// @warning Will only compile with floating point value types
315  template<typename RealT, template<typename...> class Vec3T>
316  inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
317 
318 }; // SampleFromVoxels<TreeOrAccT, 1, false>
319 
320 /// @brief Template specialization with caching of stencil values
321 template<typename TreeOrAccT>
322 class SampleFromVoxels<TreeOrAccT, 1, true> : public TrilinearSampler<TreeOrAccT>
323 {
324  using BaseT = TrilinearSampler<TreeOrAccT>;
325  using ValueT = typename TreeOrAccT::ValueType;
326  using CoordT = typename TreeOrAccT::CoordType;
327 
328  mutable CoordT mPos;
329  mutable ValueT mVal[2][2][2];
330 
331  template<typename RealT, template<typename...> class Vec3T>
332  __hostdev__ void cache(Vec3T<RealT>& xyz) const;
333 public:
334 
335  /// @brief Construction from a Tree or ReadAccessor
336  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc), mPos(CoordT::max()){}
337 
338  /// @note xyz is in index space space
339  template<typename RealT, template<typename...> class Vec3T>
340  inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
341 
342  // @note ijk is in index space space
343  __hostdev__ ValueT operator()(const CoordT &ijk) const;
344 
345  /// @brief Return the gradient in index space.
346  ///
347  /// @warning Will only compile with floating point value types
348  template<typename RealT, template<typename...> class Vec3T>
349  inline __hostdev__ Vec3T<ValueT> gradient(Vec3T<RealT> xyz) const;
350 
351  /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
352  ///
353  /// @warning Will only compile with floating point value types
354  template<typename RealT, template<typename...> class Vec3T>
355  inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
356 
357  /// @brief Return true if the cached tri-linear stencil has a zero crossing.
358  ///
359  /// @warning Will only compile with floating point value types
360  __hostdev__ bool zeroCrossing() const { return BaseT::zeroCrossing(mVal); }
361 
362 }; // SampleFromVoxels<TreeOrAccT, 1, true>
363 
364 template<typename TreeOrAccT>
365 template<typename RealT, template<typename...> class Vec3T>
366 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, true>::operator()(Vec3T<RealT> xyz) const
367 {
368  this->cache(xyz);
369  return BaseT::sample(xyz, mVal);
370 }
371 
372 template<typename TreeOrAccT>
373 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, true>::operator()(const CoordT &ijk) const
374 {
375  return ijk == mPos ? mVal[0][0][0] : BaseT::mAcc.getValue(ijk);
376 }
377 
378 template<typename TreeOrAccT>
379 template<typename RealT, template<typename...> class Vec3T>
380 __hostdev__ Vec3T<typename TreeOrAccT::ValueType> SampleFromVoxels<TreeOrAccT, 1, true>::gradient(Vec3T<RealT> xyz) const
381 {
382  this->cache(xyz);
383  return BaseT::gradient(xyz, mVal);
384 }
385 
386 template<typename TreeOrAccT>
387 template<typename RealT, template<typename...> class Vec3T>
388 __hostdev__ bool SampleFromVoxels<TreeOrAccT, 1, true>::zeroCrossing(Vec3T<RealT> xyz) const
389 {
390  this->cache(xyz);
391  return BaseT::zeroCrossing(mVal);
392 }
393 
394 template<typename TreeOrAccT>
395 template<typename RealT, template<typename...> class Vec3T>
396 __hostdev__ void SampleFromVoxels<TreeOrAccT, 1, true>::cache(Vec3T<RealT>& xyz) const
397 {
398  CoordT ijk = Floor<CoordT>(xyz);
399  if (ijk != mPos) {
400  mPos = ijk;
401  BaseT::stencil(ijk, mVal);
402  }
403 }
404 
405 #if 0
406 
407 template<typename TreeOrAccT>
408 template<typename RealT, template<typename...> class Vec3T>
409 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, false>::operator()(Vec3T<RealT> xyz) const
410 {
411  ValueT val[2][2][2];
412  CoordT ijk = Floor<CoordT>(xyz);
413  BaseT::stencil(ijk, val);
414  return BaseT::sample(xyz, val);
415 }
416 
417 #else
418 
419 template<typename TreeOrAccT>
420 template<typename RealT, template<typename...> class Vec3T>
421 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, false>::operator()(Vec3T<RealT> xyz) const
422 {
423  auto lerp = [](ValueT a, ValueT b, RealT w) { return a + ValueT(w) * (b - a); };
424 
425  CoordT coord = Floor<CoordT>(xyz);
426 
427  ValueT vx, vx1, vy, vy1, vz, vz1;
428 
429  vz = BaseT::mAcc.getValue(coord);
430  coord[2] += 1;
431  vz1 = BaseT::mAcc.getValue(coord);
432  vy = lerp(vz, vz1, xyz[2]);
433 
434  coord[1] += 1;
435 
436  vz1 = BaseT::mAcc.getValue(coord);
437  coord[2] -= 1;
438  vz = BaseT::mAcc.getValue(coord);
439  vy1 = lerp(vz, vz1, xyz[2]);
440 
441  vx = lerp(vy, vy1, xyz[1]);
442 
443  coord[0] += 1;
444 
445  vz = BaseT::mAcc.getValue(coord);
446  coord[2] += 1;
447  vz1 = BaseT::mAcc.getValue(coord);
448  vy1 = lerp(vz, vz1, xyz[2]);
449 
450  coord[1] -= 1;
451 
452  vz1 = BaseT::mAcc.getValue(coord);
453  coord[2] -= 1;
454  vz = BaseT::mAcc.getValue(coord);
455  vy = lerp(vz, vz1, xyz[2]);
456 
457  vx1 = lerp(vy, vy1, xyz[1]);
458 
459  return lerp(vx, vx1, xyz[0]);
460 }
461 #endif
462 
463 
464 template<typename TreeOrAccT>
465 template<typename RealT, template<typename...> class Vec3T>
466 __hostdev__ inline Vec3T<typename TreeOrAccT::ValueType> SampleFromVoxels<TreeOrAccT, 1, false>::gradient(Vec3T<RealT> xyz) const
467 {
468  ValueT val[2][2][2];
469  CoordT ijk = Floor<CoordT>(xyz);
470  BaseT::stencil(ijk, val);
471  return BaseT::gradient(xyz, val);
472 }
473 
474 template<typename TreeOrAccT>
475 template<typename RealT, template<typename...> class Vec3T>
476 __hostdev__ bool SampleFromVoxels<TreeOrAccT, 1, false>::zeroCrossing(Vec3T<RealT> xyz) const
477 {
478  ValueT val[2][2][2];
479  CoordT ijk = Floor<CoordT>(xyz);
480  BaseT::stencil(ijk, val);
481  return BaseT::zeroCrossing(val);
482 }
483 
484 // ------------------------------> TriquadraticSampler <--------------------------------------
485 
486 /// @brief Tri-quadratic sampler, i.e. second order, interpolator
487 template<typename TreeOrAccT>
488 class TriquadraticSampler
489 {
490 protected:
491  const TreeOrAccT& mAcc;
492 
493 public:
494  using ValueT = typename TreeOrAccT::ValueType;
495  using CoordT = typename TreeOrAccT::CoordType;
496  static const int ORDER = 1;
497 
498  /// @brief Protected constructor from a Tree or ReadAccessor
499  __hostdev__ TriquadraticSampler(const TreeOrAccT& acc) : mAcc(acc) {}
500 
501  __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
502 
503  /// @brief Extract the stencil of 27 values
504  inline __hostdev__ void stencil(const CoordT &ijk, ValueT (&v)[3][3][3]) const;
505 
506  template<typename RealT, template<typename...> class Vec3T>
507  static inline __hostdev__ ValueT sample(const Vec3T<RealT> &uvw, const ValueT (&v)[3][3][3]);
508 
509  static inline __hostdev__ bool zeroCrossing(const ValueT (&v)[3][3][3]);
510 }; // TriquadraticSamplerBase
511 
512 template<typename TreeOrAccT>
513 __hostdev__ void TriquadraticSampler<TreeOrAccT>::stencil(const CoordT &ijk, ValueT (&v)[3][3][3]) const
514 {
515  CoordT p(ijk[0] - 1, 0, 0);
516  for (int dx = 0; dx < 3; ++dx, ++p[0]) {
517  p[1] = ijk[1] - 1;
518  for (int dy = 0; dy < 3; ++dy, ++p[1]) {
519  p[2] = ijk[2] - 1;
520  for (int dz = 0; dz < 3; ++dz, ++p[2]) {
521  v[dx][dy][dz] = mAcc.getValue(p);// extract the stencil of 27 values
522  }
523  }
524  }
525 }
526 
527 template<typename TreeOrAccT>
528 template<typename RealT, template<typename...> class Vec3T>
529 __hostdev__ typename TreeOrAccT::ValueType TriquadraticSampler<TreeOrAccT>::sample(const Vec3T<RealT> &uvw, const ValueT (&v)[3][3][3])
530 {
531  auto kernel = [](const ValueT* value, double weight)->ValueT {
532  return weight * (weight * (0.5f * (value[0] + value[2]) - value[1]) +
533  0.5f * (value[2] - value[0])) + value[1];
534  };
535 
536  ValueT vx[3];
537  for (int dx = 0; dx < 3; ++dx) {
538  ValueT vy[3];
539  for (int dy = 0; dy < 3; ++dy) {
540  vy[dy] = kernel(&v[dx][dy][0], uvw[2]);
541  }//loop over y
542  vx[dx] = kernel(vy, uvw[1]);
543  }//loop over x
544  return kernel(vx, uvw[0]);
545 }
546 
547 template<typename TreeOrAccT>
548 __hostdev__ bool TriquadraticSampler<TreeOrAccT>::zeroCrossing(const ValueT (&v)[3][3][3])
549 {
550  static_assert(is_floating_point<ValueT>::value, "TrilinearSampler::zeroCrossing requires a floating-point type");
551  const bool less = v[0][0][0] < ValueT(0);
552  for (int dx = 0; dx < 3; ++dx) {
553  for (int dy = 0; dy < 3; ++dy) {
554  for (int dz = 0; dz < 3; ++dz) {
555  if (less ^ (v[dx][dy][dz] < ValueT(0))) return true;
556  }
557  }
558  }
559  return false;
560 }
561 
562 /// @brief Template specialization that does not use caching of stencil points
563 template<typename TreeOrAccT>
564 class SampleFromVoxels<TreeOrAccT, 2, false> : public TriquadraticSampler<TreeOrAccT>
565 {
566  using BaseT = TriquadraticSampler<TreeOrAccT>;
567  using ValueT = typename TreeOrAccT::ValueType;
568  using CoordT = typename TreeOrAccT::CoordType;
569 public:
570 
571  /// @brief Construction from a Tree or ReadAccessor
572  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc) {}
573 
574  /// @note xyz is in index space space
575  template<typename RealT, template<typename...> class Vec3T>
576  inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
577 
578  __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
579 
580  /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
581  ///
582  /// @warning Will only compile with floating point value types
583  template<typename RealT, template<typename...> class Vec3T>
584  inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
585 
586 }; // SampleFromVoxels<TreeOrAccT, 2, false>
587 
588 /// @brief Template specialization with caching of stencil values
589 template<typename TreeOrAccT>
590 class SampleFromVoxels<TreeOrAccT, 2, true> : public TriquadraticSampler<TreeOrAccT>
591 {
592  using BaseT = TriquadraticSampler<TreeOrAccT>;
593  using ValueT = typename TreeOrAccT::ValueType;
594  using CoordT = typename TreeOrAccT::CoordType;
595 
596  mutable CoordT mPos;
597  mutable ValueT mVal[3][3][3];
598 
599  template<typename RealT, template<typename...> class Vec3T>
600  __hostdev__ void cache(Vec3T<RealT>& xyz) const;
601 public:
602 
603  /// @brief Construction from a Tree or ReadAccessor
604  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc), mPos(CoordT::max()){}
605 
606  /// @note xyz is in index space space
607  template<typename RealT, template<typename...> class Vec3T>
608  inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
609 
610  inline __hostdev__ ValueT operator()(const CoordT &ijk) const;
611 
612  /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
613  ///
614  /// @warning Will only compile with floating point value types
615  template<typename RealT, template<typename...> class Vec3T>
616  inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
617 
618  /// @brief Return true if the cached tri-linear stencil has a zero crossing.
619  ///
620  /// @warning Will only compile with floating point value types
621  __hostdev__ bool zeroCrossing() const { return BaseT::zeroCrossing(mVal); }
622 
623 }; // SampleFromVoxels<TreeOrAccT, 2, true>
624 
625 template<typename TreeOrAccT>
626 template<typename RealT, template<typename...> class Vec3T>
627 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 2, true>::operator()(Vec3T<RealT> xyz) const
628 {
629  this->cache(xyz);
630  return BaseT::sample(xyz, mVal);
631 }
632 
633 template<typename TreeOrAccT>
634 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 2, true>::operator()(const CoordT &ijk) const
635 {
636  return ijk == mPos ? mVal[1][1][1] : BaseT::mAcc.getValue(ijk);
637 }
638 
639 template<typename TreeOrAccT>
640 template<typename RealT, template<typename...> class Vec3T>
641 __hostdev__ bool SampleFromVoxels<TreeOrAccT, 2, true>::zeroCrossing(Vec3T<RealT> xyz) const
642 {
643  this->cache(xyz);
644  return BaseT::zeroCrossing(mVal);
645 }
646 
647 template<typename TreeOrAccT>
648 template<typename RealT, template<typename...> class Vec3T>
649 __hostdev__ void SampleFromVoxels<TreeOrAccT, 2, true>::cache(Vec3T<RealT>& xyz) const
650 {
651  CoordT ijk = Floor<CoordT>(xyz);
652  if (ijk != mPos) {
653  mPos = ijk;
654  BaseT::stencil(ijk, mVal);
655  }
656 }
657 
658 template<typename TreeOrAccT>
659 template<typename RealT, template<typename...> class Vec3T>
660 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 2, false>::operator()(Vec3T<RealT> xyz) const
661 {
662  ValueT val[3][3][3];
663  CoordT ijk = Floor<CoordT>(xyz);
664  BaseT::stencil(ijk, val);
665  return BaseT::sample(xyz, val);
666 }
667 
668 template<typename TreeOrAccT>
669 template<typename RealT, template<typename...> class Vec3T>
670 __hostdev__ bool SampleFromVoxels<TreeOrAccT, 2, false>::zeroCrossing(Vec3T<RealT> xyz) const
671 {
672  ValueT val[3][3][3];
673  CoordT ijk = Floor<CoordT>(xyz);
674  BaseT::stencil(ijk, val);
675  return BaseT::zeroCrossing(val);
676 }
677 
678 // ------------------------------> TricubicSampler <--------------------------------------
679 
680 /// @brief Tri-cubic sampler, i.e. third order, interpolator.
681 ///
682 /// @details See the following paper for implementation details:
683 /// Lekien, F. and Marsden, J.: Tricubic interpolation in three dimensions.
684 /// In: International Journal for Numerical Methods
685 /// in Engineering (2005), No. 63, p. 455-471
686 
687 template<typename TreeOrAccT>
688 class TricubicSampler
689 {
690 protected:
691  using ValueT = typename TreeOrAccT::ValueType;
692  using CoordT = typename TreeOrAccT::CoordType;
693 
694  const TreeOrAccT& mAcc;
695 
696 public:
697  /// @brief Construction from a Tree or ReadAccessor
698  __hostdev__ TricubicSampler(const TreeOrAccT& acc)
699  : mAcc(acc)
700  {
701  }
702 
703  __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
704 
705  /// @brief Extract the stencil of 8 values
706  inline __hostdev__ void stencil(const CoordT& ijk, ValueT (&c)[64]) const;
707 
708  template<typename RealT, template<typename...> class Vec3T>
709  static inline __hostdev__ ValueT sample(const Vec3T<RealT> &uvw, const ValueT (&c)[64]);
710 }; // TricubicSampler
711 
712 template<typename TreeOrAccT>
713 __hostdev__ void TricubicSampler<TreeOrAccT>::stencil(const CoordT& ijk, ValueT (&C)[64]) const
714 {
715  auto fetch = [&](int i, int j, int k) -> ValueT& { return C[((i + 1) << 4) + ((j + 1) << 2) + k + 1]; };
716 
717  // fetch 64 point stencil values
718  for (int i = -1; i < 3; ++i) {
719  for (int j = -1; j < 3; ++j) {
720  fetch(i, j, -1) = mAcc.getValue(ijk + CoordT(i, j, -1));
721  fetch(i, j, 0) = mAcc.getValue(ijk + CoordT(i, j, 0));
722  fetch(i, j, 1) = mAcc.getValue(ijk + CoordT(i, j, 1));
723  fetch(i, j, 2) = mAcc.getValue(ijk + CoordT(i, j, 2));
724  }
725  }
726  const ValueT half(0.5), quarter(0.25), eighth(0.125);
727  const ValueT X[64] = {// values of f(x,y,z) at the 8 corners (each from 1 stencil value).
728  fetch(0, 0, 0),
729  fetch(1, 0, 0),
730  fetch(0, 1, 0),
731  fetch(1, 1, 0),
732  fetch(0, 0, 1),
733  fetch(1, 0, 1),
734  fetch(0, 1, 1),
735  fetch(1, 1, 1),
736  // values of df/dx at the 8 corners (each from 2 stencil values).
737  half * (fetch(1, 0, 0) - fetch(-1, 0, 0)),
738  half * (fetch(2, 0, 0) - fetch(0, 0, 0)),
739  half * (fetch(1, 1, 0) - fetch(-1, 1, 0)),
740  half * (fetch(2, 1, 0) - fetch(0, 1, 0)),
741  half * (fetch(1, 0, 1) - fetch(-1, 0, 1)),
742  half * (fetch(2, 0, 1) - fetch(0, 0, 1)),
743  half * (fetch(1, 1, 1) - fetch(-1, 1, 1)),
744  half * (fetch(2, 1, 1) - fetch(0, 1, 1)),
745  // values of df/dy at the 8 corners (each from 2 stencil values).
746  half * (fetch(0, 1, 0) - fetch(0, -1, 0)),
747  half * (fetch(1, 1, 0) - fetch(1, -1, 0)),
748  half * (fetch(0, 2, 0) - fetch(0, 0, 0)),
749  half * (fetch(1, 2, 0) - fetch(1, 0, 0)),
750  half * (fetch(0, 1, 1) - fetch(0, -1, 1)),
751  half * (fetch(1, 1, 1) - fetch(1, -1, 1)),
752  half * (fetch(0, 2, 1) - fetch(0, 0, 1)),
753  half * (fetch(1, 2, 1) - fetch(1, 0, 1)),
754  // values of df/dz at the 8 corners (each from 2 stencil values).
755  half * (fetch(0, 0, 1) - fetch(0, 0, -1)),
756  half * (fetch(1, 0, 1) - fetch(1, 0, -1)),
757  half * (fetch(0, 1, 1) - fetch(0, 1, -1)),
758  half * (fetch(1, 1, 1) - fetch(1, 1, -1)),
759  half * (fetch(0, 0, 2) - fetch(0, 0, 0)),
760  half * (fetch(1, 0, 2) - fetch(1, 0, 0)),
761  half * (fetch(0, 1, 2) - fetch(0, 1, 0)),
762  half * (fetch(1, 1, 2) - fetch(1, 1, 0)),
763  // values of d2f/dxdy at the 8 corners (each from 4 stencil values).
764  quarter * (fetch(1, 1, 0) - fetch(-1, 1, 0) - fetch(1, -1, 0) + fetch(-1, -1, 0)),
765  quarter * (fetch(2, 1, 0) - fetch(0, 1, 0) - fetch(2, -1, 0) + fetch(0, -1, 0)),
766  quarter * (fetch(1, 2, 0) - fetch(-1, 2, 0) - fetch(1, 0, 0) + fetch(-1, 0, 0)),
767  quarter * (fetch(2, 2, 0) - fetch(0, 2, 0) - fetch(2, 0, 0) + fetch(0, 0, 0)),
768  quarter * (fetch(1, 1, 1) - fetch(-1, 1, 1) - fetch(1, -1, 1) + fetch(-1, -1, 1)),
769  quarter * (fetch(2, 1, 1) - fetch(0, 1, 1) - fetch(2, -1, 1) + fetch(0, -1, 1)),
770  quarter * (fetch(1, 2, 1) - fetch(-1, 2, 1) - fetch(1, 0, 1) + fetch(-1, 0, 1)),
771  quarter * (fetch(2, 2, 1) - fetch(0, 2, 1) - fetch(2, 0, 1) + fetch(0, 0, 1)),
772  // values of d2f/dxdz at the 8 corners (each from 4 stencil values).
773  quarter * (fetch(1, 0, 1) - fetch(-1, 0, 1) - fetch(1, 0, -1) + fetch(-1, 0, -1)),
774  quarter * (fetch(2, 0, 1) - fetch(0, 0, 1) - fetch(2, 0, -1) + fetch(0, 0, -1)),
775  quarter * (fetch(1, 1, 1) - fetch(-1, 1, 1) - fetch(1, 1, -1) + fetch(-1, 1, -1)),
776  quarter * (fetch(2, 1, 1) - fetch(0, 1, 1) - fetch(2, 1, -1) + fetch(0, 1, -1)),
777  quarter * (fetch(1, 0, 2) - fetch(-1, 0, 2) - fetch(1, 0, 0) + fetch(-1, 0, 0)),
778  quarter * (fetch(2, 0, 2) - fetch(0, 0, 2) - fetch(2, 0, 0) + fetch(0, 0, 0)),
779  quarter * (fetch(1, 1, 2) - fetch(-1, 1, 2) - fetch(1, 1, 0) + fetch(-1, 1, 0)),
780  quarter * (fetch(2, 1, 2) - fetch(0, 1, 2) - fetch(2, 1, 0) + fetch(0, 1, 0)),
781  // values of d2f/dydz at the 8 corners (each from 4 stencil values).
782  quarter * (fetch(0, 1, 1) - fetch(0, -1, 1) - fetch(0, 1, -1) + fetch(0, -1, -1)),
783  quarter * (fetch(1, 1, 1) - fetch(1, -1, 1) - fetch(1, 1, -1) + fetch(1, -1, -1)),
784  quarter * (fetch(0, 2, 1) - fetch(0, 0, 1) - fetch(0, 2, -1) + fetch(0, 0, -1)),
785  quarter * (fetch(1, 2, 1) - fetch(1, 0, 1) - fetch(1, 2, -1) + fetch(1, 0, -1)),
786  quarter * (fetch(0, 1, 2) - fetch(0, -1, 2) - fetch(0, 1, 0) + fetch(0, -1, 0)),
787  quarter * (fetch(1, 1, 2) - fetch(1, -1, 2) - fetch(1, 1, 0) + fetch(1, -1, 0)),
788  quarter * (fetch(0, 2, 2) - fetch(0, 0, 2) - fetch(0, 2, 0) + fetch(0, 0, 0)),
789  quarter * (fetch(1, 2, 2) - fetch(1, 0, 2) - fetch(1, 2, 0) + fetch(1, 0, 0)),
790  // values of d3f/dxdydz at the 8 corners (each from 8 stencil values).
791  eighth * (fetch(1, 1, 1) - fetch(-1, 1, 1) - fetch(1, -1, 1) + fetch(-1, -1, 1) - fetch(1, 1, -1) + fetch(-1, 1, -1) + fetch(1, -1, -1) - fetch(-1, -1, -1)),
792  eighth * (fetch(2, 1, 1) - fetch(0, 1, 1) - fetch(2, -1, 1) + fetch(0, -1, 1) - fetch(2, 1, -1) + fetch(0, 1, -1) + fetch(2, -1, -1) - fetch(0, -1, -1)),
793  eighth * (fetch(1, 2, 1) - fetch(-1, 2, 1) - fetch(1, 0, 1) + fetch(-1, 0, 1) - fetch(1, 2, -1) + fetch(-1, 2, -1) + fetch(1, 0, -1) - fetch(-1, 0, -1)),
794  eighth * (fetch(2, 2, 1) - fetch(0, 2, 1) - fetch(2, 0, 1) + fetch(0, 0, 1) - fetch(2, 2, -1) + fetch(0, 2, -1) + fetch(2, 0, -1) - fetch(0, 0, -1)),
795  eighth * (fetch(1, 1, 2) - fetch(-1, 1, 2) - fetch(1, -1, 2) + fetch(-1, -1, 2) - fetch(1, 1, 0) + fetch(-1, 1, 0) + fetch(1, -1, 0) - fetch(-1, -1, 0)),
796  eighth * (fetch(2, 1, 2) - fetch(0, 1, 2) - fetch(2, -1, 2) + fetch(0, -1, 2) - fetch(2, 1, 0) + fetch(0, 1, 0) + fetch(2, -1, 0) - fetch(0, -1, 0)),
797  eighth * (fetch(1, 2, 2) - fetch(-1, 2, 2) - fetch(1, 0, 2) + fetch(-1, 0, 2) - fetch(1, 2, 0) + fetch(-1, 2, 0) + fetch(1, 0, 0) - fetch(-1, 0, 0)),
798  eighth * (fetch(2, 2, 2) - fetch(0, 2, 2) - fetch(2, 0, 2) + fetch(0, 0, 2) - fetch(2, 2, 0) + fetch(0, 2, 0) + fetch(2, 0, 0) - fetch(0, 0, 0))};
799 
800  // 4Kb of static table (int8_t has a range of -127 -> 127 which suffices)
801  static const int8_t A[64][64] = {
802  {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
803  {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
804  {-3, 3, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
805  {2, -2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
806  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
807  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
808  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
809  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
810  {-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
811  {0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
812  {9, -9, -9, 9, 0, 0, 0, 0, 6, 3, -6, -3, 0, 0, 0, 0, 6, -6, 3, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
813  {-6, 6, 6, -6, 0, 0, 0, 0, -3, -3, 3, 3, 0, 0, 0, 0, -4, 4, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -2, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
814  {2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
815  {0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
816  {-6, 6, 6, -6, 0, 0, 0, 0, -4, -2, 4, 2, 0, 0, 0, 0, -3, 3, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
817  {4, -4, -4, 4, 0, 0, 0, 0, 2, 2, -2, -2, 0, 0, 0, 0, 2, -2, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
818  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
819  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
820  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
821  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
822  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
823  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
824  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0},
825  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0},
826  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
827  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0},
828  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -9, -9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, -6, -3, 0, 0, 0, 0, 6, -6, 3, -3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0},
829  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, -3, 3, 3, 0, 0, 0, 0, -4, 4, -2, 2, 0, 0, 0, 0, -2, -2, -1, -1, 0, 0, 0, 0},
830  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
831  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0},
832  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, -2, 4, 2, 0, 0, 0, 0, -3, 3, -3, 3, 0, 0, 0, 0, -2, -1, -2, -1, 0, 0, 0, 0},
833  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, -4, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, -2, -2, 0, 0, 0, 0, 2, -2, 2, -2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0},
834  {-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
835  {0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
836  {9, -9, 0, 0, -9, 9, 0, 0, 6, 3, 0, 0, -6, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, -6, 0, 0, 3, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
837  {-6, 6, 0, 0, 6, -6, 0, 0, -3, -3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 4, 0, 0, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -2, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
838  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
839  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, -1, 0, 0, 0},
840  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -9, 0, 0, -9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0, -6, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, -6, 0, 0, 3, -3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0},
841  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 0, 0, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, -3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 4, 0, 0, -2, 2, 0, 0, -2, -2, 0, 0, -1, -1, 0, 0},
842  {9, 0, -9, 0, -9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0, -6, 0, -3, 0, 6, 0, -6, 0, 3, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
843  {0, 0, 0, 0, 0, 0, 0, 0, 9, 0, -9, 0, -9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0, -6, 0, -3, 0, 6, 0, -6, 0, 3, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0},
844  {-27, 27, 27, -27, 27, -27, -27, 27, -18, -9, 18, 9, 18, 9, -18, -9, -18, 18, -9, 9, 18, -18, 9, -9, -18, 18, 18, -18, -9, 9, 9, -9, -12, -6, -6, -3, 12, 6, 6, 3, -12, -6, 12, 6, -6, -3, 6, 3, -12, 12, -6, 6, -6, 6, -3, 3, -8, -4, -4, -2, -4, -2, -2, -1},
845  {18, -18, -18, 18, -18, 18, 18, -18, 9, 9, -9, -9, -9, -9, 9, 9, 12, -12, 6, -6, -12, 12, -6, 6, 12, -12, -12, 12, 6, -6, -6, 6, 6, 6, 3, 3, -6, -6, -3, -3, 6, 6, -6, -6, 3, 3, -3, -3, 8, -8, 4, -4, 4, -4, 2, -2, 4, 4, 2, 2, 2, 2, 1, 1},
846  {-6, 0, 6, 0, 6, 0, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, -3, 0, 3, 0, 3, 0, -4, 0, 4, 0, -2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
847  {0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 6, 0, 6, 0, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, -3, 0, 3, 0, 3, 0, -4, 0, 4, 0, -2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, -1, 0, -1, 0},
848  {18, -18, -18, 18, -18, 18, 18, -18, 12, 6, -12, -6, -12, -6, 12, 6, 9, -9, 9, -9, -9, 9, -9, 9, 12, -12, -12, 12, 6, -6, -6, 6, 6, 3, 6, 3, -6, -3, -6, -3, 8, 4, -8, -4, 4, 2, -4, -2, 6, -6, 6, -6, 3, -3, 3, -3, 4, 2, 4, 2, 2, 1, 2, 1},
849  {-12, 12, 12, -12, 12, -12, -12, 12, -6, -6, 6, 6, 6, 6, -6, -6, -6, 6, -6, 6, 6, -6, 6, -6, -8, 8, 8, -8, -4, 4, 4, -4, -3, -3, -3, -3, 3, 3, 3, 3, -4, -4, 4, 4, -2, -2, 2, 2, -4, 4, -4, 4, -2, 2, -2, 2, -2, -2, -2, -2, -1, -1, -1, -1},
850  {2, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
851  {0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
852  {-6, 6, 0, 0, 6, -6, 0, 0, -4, -2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
853  {4, -4, 0, 0, -4, 4, 0, 0, 2, 2, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
854  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
855  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
856  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 0, 0, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, -2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0, -2, -1, 0, 0},
857  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, -4, 0, 0, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0},
858  {-6, 0, 6, 0, 6, 0, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, -2, 0, 4, 0, 2, 0, -3, 0, 3, 0, -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
859  {0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 6, 0, 6, 0, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, -2, 0, 4, 0, 2, 0, -3, 0, 3, 0, -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, -2, 0, -1, 0},
860  {18, -18, -18, 18, -18, 18, 18, -18, 12, 6, -12, -6, -12, -6, 12, 6, 12, -12, 6, -6, -12, 12, -6, 6, 9, -9, -9, 9, 9, -9, -9, 9, 8, 4, 4, 2, -8, -4, -4, -2, 6, 3, -6, -3, 6, 3, -6, -3, 6, -6, 3, -3, 6, -6, 3, -3, 4, 2, 2, 1, 4, 2, 2, 1},
861  {-12, 12, 12, -12, 12, -12, -12, 12, -6, -6, 6, 6, 6, 6, -6, -6, -8, 8, -4, 4, 8, -8, 4, -4, -6, 6, 6, -6, -6, 6, 6, -6, -4, -4, -2, -2, 4, 4, 2, 2, -3, -3, 3, 3, -3, -3, 3, 3, -4, 4, -2, 2, -4, 4, -2, 2, -2, -2, -1, -1, -2, -2, -1, -1},
862  {4, 0, -4, 0, -4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0, -2, 0, -2, 0, 2, 0, -2, 0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
863  {0, 0, 0, 0, 0, 0, 0, 0, 4, 0, -4, 0, -4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0, -2, 0, -2, 0, 2, 0, -2, 0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0},
864  {-12, 12, 12, -12, 12, -12, -12, 12, -8, -4, 8, 4, 8, 4, -8, -4, -6, 6, -6, 6, 6, -6, 6, -6, -6, 6, 6, -6, -6, 6, 6, -6, -4, -2, -4, -2, 4, 2, 4, 2, -4, -2, 4, 2, -4, -2, 4, 2, -3, 3, -3, 3, -3, 3, -3, 3, -2, -1, -2, -1, -2, -1, -2, -1},
865  {8, -8, -8, 8, -8, 8, 8, -8, 4, 4, -4, -4, -4, -4, 4, 4, 4, -4, 4, -4, -4, 4, -4, 4, 4, -4, -4, 4, 4, -4, -4, 4, 2, 2, 2, 2, -2, -2, -2, -2, 2, 2, -2, -2, 2, 2, -2, -2, 2, -2, 2, -2, 2, -2, 2, -2, 1, 1, 1, 1, 1, 1, 1, 1}};
866 
867  for (int i = 0; i < 64; ++i) { // C = A * X
868  C[i] = ValueT(0);
869 #if 0
870  for (int j = 0; j < 64; j += 4) {
871  C[i] = fma(A[i][j], X[j], fma(A[i][j+1], X[j+1], fma(A[i][j+2], X[j+2], fma(A[i][j+3], X[j+3], C[i]))));
872  }
873 #else
874  for (int j = 0; j < 64; j += 4) {
875  C[i] += A[i][j] * X[j] + A[i][j + 1] * X[j + 1] + A[i][j + 2] * X[j + 2] + A[i][j + 3] * X[j + 3];
876  }
877 #endif
878  }
879 }
880 
881 template<typename TreeOrAccT>
882 template<typename RealT, template<typename...> class Vec3T>
883 __hostdev__ typename TreeOrAccT::ValueType TricubicSampler<TreeOrAccT>::sample(const Vec3T<RealT> &xyz, const ValueT (&C)[64])
884 {
885  ValueT zPow(1), sum(0);
886  for (int k = 0, n = 0; k < 4; ++k) {
887  ValueT yPow(1);
888  for (int j = 0; j < 4; ++j, n += 4) {
889 #if 0
890  sum = fma( yPow, zPow * fma(xyz[0], fma(xyz[0], fma(xyz[0], C[n + 3], C[n + 2]), C[n + 1]), C[n]), sum);
891 #else
892  sum += yPow * zPow * (C[n] + xyz[0] * (C[n + 1] + xyz[0] * (C[n + 2] + xyz[0] * C[n + 3])));
893 #endif
894  yPow *= xyz[1];
895  }
896  zPow *= xyz[2];
897  }
898  return sum;
899 }
900 
901 template<typename TreeOrAccT>
902 class SampleFromVoxels<TreeOrAccT, 3, true> : public TricubicSampler<TreeOrAccT>
903 {
904  using BaseT = TricubicSampler<TreeOrAccT>;
905  using ValueT = typename TreeOrAccT::ValueType;
906  using CoordT = typename TreeOrAccT::CoordType;
907 
908  mutable CoordT mPos;
909  mutable ValueT mC[64];
910 
911  template<typename RealT, template<typename...> class Vec3T>
912  __hostdev__ void cache(Vec3T<RealT>& xyz) const;
913 
914 public:
915  /// @brief Construction from a Tree or ReadAccessor
916  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc)
917  : BaseT(acc)
918  {
919  }
920 
921  /// @note xyz is in index space space
922  template<typename RealT, template<typename...> class Vec3T>
923  inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
924 
925  // @brief Return value at the coordinate @a ijk in index space space
926  __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
927 
928 }; // SampleFromVoxels<TreeOrAccT, 3, true>
929 
930 template<typename TreeOrAccT>
931 template<typename RealT, template<typename...> class Vec3T>
932 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 3, true>::operator()(Vec3T<RealT> xyz) const
933 {
934  this->cache(xyz);
935  return BaseT::sample(xyz, mC);
936 }
937 
938 template<typename TreeOrAccT>
939 template<typename RealT, template<typename...> class Vec3T>
940 __hostdev__ void SampleFromVoxels<TreeOrAccT, 3, true>::cache(Vec3T<RealT>& xyz) const
941 {
942  CoordT ijk = Floor<CoordT>(xyz);
943  if (ijk != mPos) {
944  mPos = ijk;
945  BaseT::stencil(ijk, mC);
946  }
947 }
948 
949 template<typename TreeOrAccT>
950 class SampleFromVoxels<TreeOrAccT, 3, false> : public TricubicSampler<TreeOrAccT>
951 {
952  using BaseT = TricubicSampler<TreeOrAccT>;
953  using ValueT = typename TreeOrAccT::ValueType;
954  using CoordT = typename TreeOrAccT::CoordType;
955 
956 public:
957  /// @brief Construction from a Tree or ReadAccessor
958  __hostdev__ SampleFromVoxels(const TreeOrAccT& acc)
959  : BaseT(acc)
960  {
961  }
962 
963  /// @note xyz is in index space space
964  template<typename RealT, template<typename...> class Vec3T>
965  inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
966 
967  __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
968 
969 }; // SampleFromVoxels<TreeOrAccT, 3, true>
970 
971 template<typename TreeOrAccT>
972 template<typename RealT, template<typename...> class Vec3T>
973 __hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 3, false>::operator()(Vec3T<RealT> xyz) const
974 {
975  ValueT C[64];
976  CoordT ijk = Floor<CoordT>(xyz);
977  BaseT::stencil(ijk, C);
978  return BaseT::sample(xyz, C);
979 }
980 
981 } // namespace nanovdb
982 
983 #endif // NANOVDB_SAMPLE_FROM_VOXELS_H_HAS_BEEN_INCLUDED
nanovdb::TrilinearSampler::accessor
__hostdev__ const TreeOrAccT & accessor() const
Definition: SampleFromVoxels.h:183
nanovdb::TrilinearSampler::ValueT
typename TreeOrAccT::ValueType ValueT
Definition: SampleFromVoxels.h:176
nanovdb::SampleFromVoxels< TreeOrAccT, 1, true >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:336
nanovdb::SampleFromVoxels< TreeOrAccT, 1, false >::operator()
__hostdev__ ValueT operator()(const CoordT &ijk) const
Definition: SampleFromVoxels.h:304
__hostdev__
#define __hostdev__
Definition: SampleFromVoxels.h:29
nanovdb::createSampler
__hostdev__ SampleFromVoxels< TreeOrAccT, Order, UseCache > createSampler(const TreeOrAccT &acc)
Factory free-function for a sampler of specific polynomial orders.
Definition: SampleFromVoxels.h:46
half
imath_half_bits_t half
if we're in a C-only context, alias the half bits type to half
Definition: half.h:266
nanovdb::SampleFromVoxels< TreeOrAccT, 2, false >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:572
nanovdb::TriquadraticSampler::ValueT
typename TreeOrAccT::ValueType ValueT
Definition: SampleFromVoxels.h:494
floor
IMATH_HOSTDEVICE constexpr int floor(T x) IMATH_NOEXCEPT
Definition: ImathFun.h:112
nanovdb::TriquadraticSampler::stencil
__hostdev__ void stencil(const CoordT &ijk, ValueT(&v)[3][3][3]) const
Extract the stencil of 27 values.
Definition: SampleFromVoxels.h:513
v
const GLdouble * v
Definition: glcorearb.h:837
value
GLsizei const GLfloat * value
Definition: glcorearb.h:824
nanovdb::TricubicSampler::stencil
__hostdev__ void stencil(const CoordT &ijk, ValueT(&c)[64]) const
Extract the stencil of 8 values.
Definition: SampleFromVoxels.h:713
nanovdb::TrilinearSampler::zeroCrossing
static __hostdev__ bool zeroCrossing(const ValueT(&v)[2][2][2])
Definition: SampleFromVoxels.h:273
lerp
virtual bool lerp(GA_AttributeOperand &d, GA_AttributeOperand &a, GA_AttributeOperand &b, GA_AttributeOperand &t) const
d = SYSlerp(a, b, t);
nanovdb::TriquadraticSampler
Tri-quadratic sampler, i.e. second order, interpolator.
Definition: SampleFromVoxels.h:488
nanovdb::TricubicSampler::CoordT
typename TreeOrAccT::CoordType CoordT
Definition: SampleFromVoxels.h:692
a
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
X
X
Definition: ImathEuler.h:183
nanovdb::TriquadraticSampler::sample
static __hostdev__ ValueT sample(const Vec3T< RealT > &uvw, const ValueT(&v)[3][3][3])
nanovdb::SampleFromVoxels< TreeOrAccT, 3, false >::operator()
__hostdev__ ValueT operator()(const CoordT &ijk) const
Definition: SampleFromVoxels.h:967
nanovdb::TricubicSampler::sample
static __hostdev__ ValueT sample(const Vec3T< RealT > &uvw, const ValueT(&c)[64])
nanovdb::SampleFromVoxels< TreeOrAccT, 2, true >::zeroCrossing
__hostdev__ bool zeroCrossing() const
Return true if the cached tri-linear stencil has a zero crossing.
Definition: SampleFromVoxels.h:621
nanovdb::TrilinearSampler::TrilinearSampler
__hostdev__ TrilinearSampler(const TreeOrAccT &acc)
Protected constructor from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:181
nanovdb::TriquadraticSampler::zeroCrossing
static __hostdev__ bool zeroCrossing(const ValueT(&v)[3][3][3])
Definition: SampleFromVoxels.h:548
nanovdb::TricubicSampler::ValueT
typename TreeOrAccT::ValueType ValueT
Definition: SampleFromVoxels.h:691
nanovdb::SampleFromVoxels< TreeOrAccT, 1, true >::zeroCrossing
__hostdev__ bool zeroCrossing() const
Return true if the cached tri-linear stencil has a zero crossing.
Definition: SampleFromVoxels.h:360
n
GLdouble n
Definition: glcorearb.h:2008
nanovdb::SampleFromVoxels< TreeOrAccT, 3, true >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:916
nanovdb::TricubicSampler
Tri-cubic sampler, i.e. third order, interpolator.
Definition: SampleFromVoxels.h:688
nanovdb::SampleFromVoxels< TreeOrAccT, 1, false >::operator()
__hostdev__ ValueT operator()(Vec3T< RealT > xyz) const
nanovdb::TrilinearSampler::stencil
__hostdev__ void stencil(CoordT &ijk, ValueT(&v)[2][2][2]) const
Extract the stencil of 8 values.
Definition: SampleFromVoxels.h:198
nanovdb::SampleFromVoxels< TreeOrAccT, 0, true >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:90
nanovdb::SampleFromVoxels< TreeOrAccT, 0, false >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:120
b
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
nanovdb::TrilinearSampler::sample
static __hostdev__ ValueT sample(const Vec3T< RealT > &uvw, const ValueT(&v)[2][2][2])
nanovdb::TrilinearSampler::gradient
static __hostdev__ Vec3T< ValueT > gradient(const Vec3T< RealT > &uvw, const ValueT(&v)[2][2][2])
j
GLint j
Definition: glad.h:2733
nanovdb::SampleFromVoxels< TreeOrAccT, 0, true >::accessor
__hostdev__ const TreeOrAccT & accessor() const
Definition: SampleFromVoxels.h:96
nanovdb::SampleFromVoxels< TreeOrAccT, 1, false >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:297
nanovdb::TriquadraticSampler::CoordT
typename TreeOrAccT::CoordType CoordT
Definition: SampleFromVoxels.h:495
nanovdb::TricubicSampler::TricubicSampler
__hostdev__ TricubicSampler(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:698
openvdb::OPENVDB_VERSION_NAME::tools::gradient
ScalarToVectorConverter< GridType >::Type::Ptr gradient(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the gradient of the given scalar grid.
Definition: GridOperators.h:1001
nanovdb::SampleFromVoxels< TreeOrAccT, 2, true >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:604
UT_Gallery::A
Definition: UT_GalleryTypes.h:40
nanovdb::Floor
__hostdev__ int32_t Floor(float x)
Definition: NanoVDB.h:1130
val
GLuint GLfloat * val
Definition: glcorearb.h:1608
ImageBufAlgo::max
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
nanovdb::SampleFromVoxels< TreeOrAccT, 3, false >::SampleFromVoxels
__hostdev__ SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:958
nanovdb::TrilinearSampler
Tri-linear sampler, i.e. first order, interpolator.
Definition: SampleFromVoxels.h:170
nanovdb::SampleFromVoxels< TreeOrAccT, 0, false >::accessor
__hostdev__ const TreeOrAccT & accessor() const
Definition: SampleFromVoxels.h:125
nanovdb::SampleFromVoxels< TreeOrAccT, 0, false >::operator()
__hostdev__ ValueT operator()(const CoordT &ijk) const
Definition: SampleFromVoxels.h:131
w
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:857
nanovdb::TricubicSampler::mAcc
const TreeOrAccT & mAcc
Definition: SampleFromVoxels.h:694
nanovdb::TriquadraticSampler::accessor
__hostdev__ const TreeOrAccT & accessor() const
Definition: SampleFromVoxels.h:501
nanovdb::SampleFromVoxels< TreeOrAccT, 3, true >::operator()
__hostdev__ ValueT operator()(const CoordT &ijk) const
Definition: SampleFromVoxels.h:926
nanovdb::TrilinearSampler::CoordT
typename TreeOrAccT::CoordType CoordT
Definition: SampleFromVoxels.h:177
stencil
GLint GLfloat GLint stencil
Definition: glcorearb.h:1278
nanovdb::TricubicSampler::accessor
__hostdev__ const TreeOrAccT & accessor() const
Definition: SampleFromVoxels.h:703
nanovdb::ValueType
bool ValueType
Definition: NanoVDB.h:5729
nanovdb::SampleFromVoxels< TreeOrAccT, 2, false >::operator()
__hostdev__ ValueT operator()(const CoordT &ijk) const
Definition: SampleFromVoxels.h:578
nanovdb::is_floating_point
C++11 implementation of std::is_floating_point.
Definition: NanoVDB.h:439
nanovdb::TriquadraticSampler::TriquadraticSampler
__hostdev__ TriquadraticSampler(const TreeOrAccT &acc)
Protected constructor from a Tree or ReadAccessor.
Definition: SampleFromVoxels.h:499