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