OpenVDB 14.0.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
172namespace detail {
173/// @brief Scalar component type of a scalar or vector value type, e.g. float for both float and Vec3<float>
174template<typename...>
175struct MakeVoid { using type = void; };// like std::void_t, which is unavailable under NVRTC
176template<typename T, typename = void>
177struct ScalarOf { using type = T; };
178template<typename T>
179struct ScalarOf<T, typename MakeVoid<typename T::ValueType>::type> { using type = typename T::ValueType; };
180
181/// @brief Explicitly convert an integral @a v to RealT, i.e. the implicit promotion in mixed
182/// integer/floating-point arithmetic. All other types are returned unchanged.
183/// @note The samplers interpolate integral grids in integer arithmetic and truncate the result,
184/// e.g. the linear weights and tricubic stencil weights become zero, so their output
185/// is generally not meaningful.
186template<typename RealT, typename T>
187__hostdev__ inline constexpr auto toReal(const T& v)
188{
189 if constexpr(std::is_integral<T>::value) {
190 return RealT(v);
191 } else {
192 return v;
193 }
194}
195
196/// @brief Preserve NanoVDB vectors' component-precision scalar multiplication.
197/// For other types, multiply before converting the result back to T.
198template<typename T, typename S>
199__hostdev__ inline T scale(const T& v, S s)
200{
201 using ScalarT = typename ScalarOf<T>::type;
202 if constexpr(util::is_specialization<T, Vec3>::value ||
203 util::is_specialization<T, Vec4>::value) {
204 return v * ScalarT(s);
205 } else {
206 return T(toReal<S>(v) * s);
207 }
208}
209}// namespace detail
210
211/// @brief Tri-linear sampler, i.e. first order, interpolator
212template<typename TreeOrAccT>
213class TrilinearSampler
214{
215protected:
216 const TreeOrAccT& mAcc;
217
218public:
219 using ValueT = typename TreeOrAccT::ValueType;
220 using CoordT = typename TreeOrAccT::CoordType;
221 static const int ORDER = 1;
222
223 /// @brief Protected constructor from a Tree or ReadAccessor
224 __hostdev__ TrilinearSampler(const TreeOrAccT& acc) : mAcc(acc) {}
225
226 __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
227
228 /// @brief Extract the stencil of 8 values
229 inline __hostdev__ void stencil(CoordT& ijk, ValueT (&v)[2][2][2]) const;
230
231 template<typename RealT, template<typename...> class Vec3T>
232 static inline __hostdev__ ValueT sample(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2]);
233
234 template<typename RealT, template<typename...> class Vec3T>
235 static inline __hostdev__ Vec3T<ValueT> gradient(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2]);
236
237 static inline __hostdev__ bool zeroCrossing(const ValueT (&v)[2][2][2]);
238}; // TrilinearSamplerBase
239
240template<typename TreeOrAccT>
241__hostdev__ void TrilinearSampler<TreeOrAccT>::stencil(CoordT& ijk, ValueT (&v)[2][2][2]) const
242{
243 v[0][0][0] = mAcc.getValue(ijk); // i, j, k
244
245 ijk[2] += 1;
246 v[0][0][1] = mAcc.getValue(ijk); // i, j, k + 1
247
248 ijk[1] += 1;
249 v[0][1][1] = mAcc.getValue(ijk); // i, j+1, k + 1
250
251 ijk[2] -= 1;
252 v[0][1][0] = mAcc.getValue(ijk); // i, j+1, k
253
254 ijk[0] += 1;
255 ijk[1] -= 1;
256 v[1][0][0] = mAcc.getValue(ijk); // i+1, j, k
257
258 ijk[2] += 1;
259 v[1][0][1] = mAcc.getValue(ijk); // i+1, j, k + 1
260
261 ijk[1] += 1;
262 v[1][1][1] = mAcc.getValue(ijk); // i+1, j+1, k + 1
263
264 ijk[2] -= 1;
265 v[1][1][0] = mAcc.getValue(ijk); // i+1, j+1, k
266}
267
268template<typename TreeOrAccT>
269template<typename RealT, template<typename...> class Vec3T>
270__hostdev__ typename TreeOrAccT::ValueType TrilinearSampler<TreeOrAccT>::sample(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2])
271{
272#if 0
273 auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b-a, a); };// = w*(b-a) + a
274 //auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b, fma(-w, a, a));};// = (1-w)*a + w*b
275#else
276 using ScalarT = typename detail::ScalarOf<ValueT>::type;
277 auto lerp = [](ValueT a, ValueT b, RealT w) { return a + ValueT(ScalarT(w)) * (b - a); };
278#endif
279 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]),
280 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]),
281 uvw[0]);
282}
283
284template<typename TreeOrAccT>
285template<typename RealT, template<typename...> class Vec3T>
286__hostdev__ Vec3T<typename TreeOrAccT::ValueType> TrilinearSampler<TreeOrAccT>::gradient(const Vec3T<RealT> &uvw, const ValueT (&v)[2][2][2])
287{
288 static_assert(util::is_floating_point<ValueT>::value, "TrilinearSampler::gradient requires a floating-point type");
289#if 0
290 auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b-a, a); };// = w*(b-a) + a
291 //auto lerp = [](ValueT a, ValueT b, ValueT w){ return fma(w, b, fma(-w, a, a));};// = (1-w)*a + w*b
292#else
293 auto lerp = [](ValueT a, ValueT b, RealT w) { return a + ValueT(w) * (b - a); };
294#endif
295
296 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]};
297
298 // Z component
299 Vec3T<ValueT> grad(0, 0, lerp(lerp(D[0], D[1], uvw[1]), lerp(D[2], D[3], uvw[1]), uvw[0]));
300
301 const ValueT w = ValueT(uvw[2]);
302 D[0] = v[0][0][0] + D[0] * w;
303 D[1] = v[0][1][0] + D[1] * w;
304 D[2] = v[1][0][0] + D[2] * w;
305 D[3] = v[1][1][0] + D[3] * w;
306
307 // X component
308 grad[0] = lerp(D[2], D[3], uvw[1]) - lerp(D[0], D[1], uvw[1]);
309
310 // Y component
311 grad[1] = lerp(D[1] - D[0], D[3] - D[2], uvw[0]);
312
313 return grad;
314}
315
316template<typename TreeOrAccT>
317__hostdev__ bool TrilinearSampler<TreeOrAccT>::zeroCrossing(const ValueT (&v)[2][2][2])
318{
319 static_assert(util::is_floating_point<ValueT>::value, "TrilinearSampler::zeroCrossing requires a floating-point type");
320 const bool less = v[0][0][0] < ValueT(0);
321 return (less ^ (v[0][0][1] < ValueT(0))) ||
322 (less ^ (v[0][1][1] < ValueT(0))) ||
323 (less ^ (v[0][1][0] < ValueT(0))) ||
324 (less ^ (v[1][0][0] < ValueT(0))) ||
325 (less ^ (v[1][0][1] < ValueT(0))) ||
326 (less ^ (v[1][1][1] < ValueT(0))) ||
327 (less ^ (v[1][1][0] < ValueT(0)));
328}
329
330/// @brief Template specialization that does not use caching of stencil points
331template<typename TreeOrAccT>
332class SampleFromVoxels<TreeOrAccT, 1, false> : public TrilinearSampler<TreeOrAccT>
333{
334 using BaseT = TrilinearSampler<TreeOrAccT>;
335 using ValueT = typename TreeOrAccT::ValueType;
336 using CoordT = typename TreeOrAccT::CoordType;
337
338public:
339
340 /// @brief Construction from a Tree or ReadAccessor
341 __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc) {}
342
343 /// @note xyz is in index space space
344 template<typename RealT, template<typename...> class Vec3T>
345 inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
346
347 /// @note ijk is in index space space
348 __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
349
350 /// @brief Return the gradient in index space.
351 ///
352 /// @warning Will only compile with floating point value types
353 template<typename RealT, template<typename...> class Vec3T>
354 inline __hostdev__ Vec3T<ValueT> gradient(Vec3T<RealT> xyz) const;
355
356 /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
357 ///
358 /// @warning Will only compile with floating point value types
359 template<typename RealT, template<typename...> class Vec3T>
360 inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
361
362}; // SampleFromVoxels<TreeOrAccT, 1, false>
363
364/// @brief Template specialization with caching of stencil values
365template<typename TreeOrAccT>
366class SampleFromVoxels<TreeOrAccT, 1, true> : public TrilinearSampler<TreeOrAccT>
367{
368 using BaseT = TrilinearSampler<TreeOrAccT>;
369 using ValueT = typename TreeOrAccT::ValueType;
370 using CoordT = typename TreeOrAccT::CoordType;
371
372 mutable CoordT mPos;
373 mutable ValueT mVal[2][2][2];
374
375 template<typename RealT, template<typename...> class Vec3T>
376 __hostdev__ void cache(Vec3T<RealT>& xyz) const;
377public:
378
379 /// @brief Construction from a Tree or ReadAccessor
380 __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc), mPos(CoordT::max()){}
381
382 /// @note xyz is in index space space
383 template<typename RealT, template<typename...> class Vec3T>
384 inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
385
386 // @note ijk is in index space space
387 __hostdev__ ValueT operator()(const CoordT &ijk) const;
388
389 /// @brief Return the gradient in index space.
390 ///
391 /// @warning Will only compile with floating point value types
392 template<typename RealT, template<typename...> class Vec3T>
393 inline __hostdev__ Vec3T<ValueT> gradient(Vec3T<RealT> xyz) const;
394
395 /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
396 ///
397 /// @warning Will only compile with floating point value types
398 template<typename RealT, template<typename...> class Vec3T>
399 inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
400
401 /// @brief Return true if the cached tri-linear stencil has a zero crossing.
402 ///
403 /// @warning Will only compile with floating point value types
404 __hostdev__ bool zeroCrossing() const { return BaseT::zeroCrossing(mVal); }
405
406}; // SampleFromVoxels<TreeOrAccT, 1, true>
407
408template<typename TreeOrAccT>
409template<typename RealT, template<typename...> class Vec3T>
410__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, true>::operator()(Vec3T<RealT> xyz) const
411{
412 this->cache(xyz);
413 return BaseT::sample(xyz, mVal);
414}
415
416template<typename TreeOrAccT>
417__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, true>::operator()(const CoordT &ijk) const
418{
419 return ijk == mPos ? mVal[0][0][0] : BaseT::mAcc.getValue(ijk);
420}
421
422template<typename TreeOrAccT>
423template<typename RealT, template<typename...> class Vec3T>
424__hostdev__ Vec3T<typename TreeOrAccT::ValueType> SampleFromVoxels<TreeOrAccT, 1, true>::gradient(Vec3T<RealT> xyz) const
425{
426 this->cache(xyz);
427 return BaseT::gradient(xyz, mVal);
428}
429
430template<typename TreeOrAccT>
431template<typename RealT, template<typename...> class Vec3T>
432__hostdev__ bool SampleFromVoxels<TreeOrAccT, 1, true>::zeroCrossing(Vec3T<RealT> xyz) const
433{
434 this->cache(xyz);
435 return BaseT::zeroCrossing(mVal);
436}
437
438template<typename TreeOrAccT>
439template<typename RealT, template<typename...> class Vec3T>
440__hostdev__ void SampleFromVoxels<TreeOrAccT, 1, true>::cache(Vec3T<RealT>& xyz) const
441{
442 CoordT ijk = Floor<CoordT>(xyz);
443 if (ijk != mPos) {
444 mPos = ijk;
445 BaseT::stencil(ijk, mVal);
446 }
447}
448
449#if 0
450
451template<typename TreeOrAccT>
452template<typename RealT, template<typename...> class Vec3T>
453__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, false>::operator()(Vec3T<RealT> xyz) const
454{
455 ValueT val[2][2][2];
456 CoordT ijk = Floor<CoordT>(xyz);
457 BaseT::stencil(ijk, val);
458 return BaseT::sample(xyz, val);
459}
460
461#else
462
463template<typename TreeOrAccT>
464template<typename RealT, template<typename...> class Vec3T>
465__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 1, false>::operator()(Vec3T<RealT> xyz) const
466{
467 using ScalarT = typename detail::ScalarOf<ValueT>::type;
468 auto lerp = [](ValueT a, ValueT b, RealT w) { return a + ValueT(ScalarT(w)) * (b - a); };
469
470 CoordT coord = Floor<CoordT>(xyz);
471
472 ValueT vx, vx1, vy, vy1, vz, vz1;
473
474 vz = BaseT::mAcc.getValue(coord);
475 coord[2] += 1;
476 vz1 = BaseT::mAcc.getValue(coord);
477 vy = lerp(vz, vz1, xyz[2]);
478
479 coord[1] += 1;
480
481 vz1 = BaseT::mAcc.getValue(coord);
482 coord[2] -= 1;
483 vz = BaseT::mAcc.getValue(coord);
484 vy1 = lerp(vz, vz1, xyz[2]);
485
486 vx = lerp(vy, vy1, xyz[1]);
487
488 coord[0] += 1;
489
490 vz = BaseT::mAcc.getValue(coord);
491 coord[2] += 1;
492 vz1 = BaseT::mAcc.getValue(coord);
493 vy1 = lerp(vz, vz1, xyz[2]);
494
495 coord[1] -= 1;
496
497 vz1 = BaseT::mAcc.getValue(coord);
498 coord[2] -= 1;
499 vz = BaseT::mAcc.getValue(coord);
500 vy = lerp(vz, vz1, xyz[2]);
501
502 vx1 = lerp(vy, vy1, xyz[1]);
503
504 return lerp(vx, vx1, xyz[0]);
505}
506#endif
507
508
509template<typename TreeOrAccT>
510template<typename RealT, template<typename...> class Vec3T>
511__hostdev__ inline Vec3T<typename TreeOrAccT::ValueType> SampleFromVoxels<TreeOrAccT, 1, false>::gradient(Vec3T<RealT> xyz) const
512{
513 ValueT val[2][2][2];
514 CoordT ijk = Floor<CoordT>(xyz);
515 BaseT::stencil(ijk, val);
516 return BaseT::gradient(xyz, val);
517}
518
519template<typename TreeOrAccT>
520template<typename RealT, template<typename...> class Vec3T>
521__hostdev__ bool SampleFromVoxels<TreeOrAccT, 1, false>::zeroCrossing(Vec3T<RealT> xyz) const
522{
523 ValueT val[2][2][2];
524 CoordT ijk = Floor<CoordT>(xyz);
525 BaseT::stencil(ijk, val);
526 return BaseT::zeroCrossing(val);
527}
528
529// ------------------------------> TriquadraticSampler <--------------------------------------
530
531/// @brief Tri-quadratic sampler, i.e. second order, interpolator
532template<typename TreeOrAccT>
533class TriquadraticSampler
534{
535protected:
536 const TreeOrAccT& mAcc;
537
538public:
539 using ValueT = typename TreeOrAccT::ValueType;
540 using CoordT = typename TreeOrAccT::CoordType;
541 static const int ORDER = 1;
542
543 /// @brief Protected constructor from a Tree or ReadAccessor
544 __hostdev__ TriquadraticSampler(const TreeOrAccT& acc) : mAcc(acc) {}
545
546 __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
547
548 /// @brief Extract the stencil of 27 values
549 inline __hostdev__ void stencil(const CoordT &ijk, ValueT (&v)[3][3][3]) const;
550
551 template<typename RealT, template<typename...> class Vec3T>
552 static inline __hostdev__ ValueT sample(const Vec3T<RealT> &uvw, const ValueT (&v)[3][3][3]);
553
554 static inline __hostdev__ bool zeroCrossing(const ValueT (&v)[3][3][3]);
555}; // TriquadraticSamplerBase
556
557template<typename TreeOrAccT>
558__hostdev__ void TriquadraticSampler<TreeOrAccT>::stencil(const CoordT &ijk, ValueT (&v)[3][3][3]) const
559{
560 CoordT p(ijk[0] - 1, 0, 0);
561 for (int dx = 0; dx < 3; ++dx, ++p[0]) {
562 p[1] = ijk[1] - 1;
563 for (int dy = 0; dy < 3; ++dy, ++p[1]) {
564 p[2] = ijk[2] - 1;
565 for (int dz = 0; dz < 3; ++dz, ++p[2]) {
566 v[dx][dy][dz] = mAcc.getValue(p);// extract the stencil of 27 values
567 }
568 }
569 }
570}
571
572template<typename TreeOrAccT>
573template<typename RealT, template<typename...> class Vec3T>
574__hostdev__ typename TreeOrAccT::ValueType TriquadraticSampler<TreeOrAccT>::sample(const Vec3T<RealT> &uvw, const ValueT (&v)[3][3][3])
575{
576 auto kernel = [](const ValueT* value, double weight)->ValueT {
577 using detail::toReal;
578 return ValueT(weight * (weight * (0.5f * toReal<float>(value[0] + value[2]) - toReal<float>(value[1])) +
579 0.5f * toReal<float>(value[2] - value[0])) + value[1]);
580 };
581
582 ValueT vx[3];
583 for (int dx = 0; dx < 3; ++dx) {
584 ValueT vy[3];
585 for (int dy = 0; dy < 3; ++dy) {
586 vy[dy] = kernel(&v[dx][dy][0], uvw[2]);
587 }//loop over y
588 vx[dx] = kernel(vy, uvw[1]);
589 }//loop over x
590 return kernel(vx, uvw[0]);
591}
592
593template<typename TreeOrAccT>
594__hostdev__ bool TriquadraticSampler<TreeOrAccT>::zeroCrossing(const ValueT (&v)[3][3][3])
595{
596 static_assert(util::is_floating_point<ValueT>::value, "TrilinearSampler::zeroCrossing requires a floating-point type");
597 const bool less = v[0][0][0] < ValueT(0);
598 for (int dx = 0; dx < 3; ++dx) {
599 for (int dy = 0; dy < 3; ++dy) {
600 for (int dz = 0; dz < 3; ++dz) {
601 if (less ^ (v[dx][dy][dz] < ValueT(0))) return true;
602 }
603 }
604 }
605 return false;
606}
607
608/// @brief Template specialization that does not use caching of stencil points
609template<typename TreeOrAccT>
610class SampleFromVoxels<TreeOrAccT, 2, false> : public TriquadraticSampler<TreeOrAccT>
611{
612 using BaseT = TriquadraticSampler<TreeOrAccT>;
613 using ValueT = typename TreeOrAccT::ValueType;
614 using CoordT = typename TreeOrAccT::CoordType;
615public:
616
617 /// @brief Construction from a Tree or ReadAccessor
618 __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc) {}
619
620 /// @note xyz is in index space space
621 template<typename RealT, template<typename...> class Vec3T>
622 inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
623
624 __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
625
626 /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
627 ///
628 /// @warning Will only compile with floating point value types
629 template<typename RealT, template<typename...> class Vec3T>
630 inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
631
632}; // SampleFromVoxels<TreeOrAccT, 2, false>
633
634/// @brief Template specialization with caching of stencil values
635template<typename TreeOrAccT>
636class SampleFromVoxels<TreeOrAccT, 2, true> : public TriquadraticSampler<TreeOrAccT>
637{
638 using BaseT = TriquadraticSampler<TreeOrAccT>;
639 using ValueT = typename TreeOrAccT::ValueType;
640 using CoordT = typename TreeOrAccT::CoordType;
641
642 mutable CoordT mPos;
643 mutable ValueT mVal[3][3][3];
644
645 template<typename RealT, template<typename...> class Vec3T>
646 __hostdev__ void cache(Vec3T<RealT>& xyz) const;
647public:
648
649 /// @brief Construction from a Tree or ReadAccessor
650 __hostdev__ SampleFromVoxels(const TreeOrAccT& acc) : BaseT(acc), mPos(CoordT::max()){}
651
652 /// @note xyz is in index space space
653 template<typename RealT, template<typename...> class Vec3T>
654 inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
655
656 inline __hostdev__ ValueT operator()(const CoordT &ijk) const;
657
658 /// @brief Return true if the tr-linear stencil has a zero crossing at the specified index position.
659 ///
660 /// @warning Will only compile with floating point value types
661 template<typename RealT, template<typename...> class Vec3T>
662 inline __hostdev__ bool zeroCrossing(Vec3T<RealT> xyz) const;
663
664 /// @brief Return true if the cached tri-linear stencil has a zero crossing.
665 ///
666 /// @warning Will only compile with floating point value types
667 __hostdev__ bool zeroCrossing() const { return BaseT::zeroCrossing(mVal); }
668
669}; // SampleFromVoxels<TreeOrAccT, 2, true>
670
671template<typename TreeOrAccT>
672template<typename RealT, template<typename...> class Vec3T>
673__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 2, true>::operator()(Vec3T<RealT> xyz) const
674{
675 this->cache(xyz);
676 return BaseT::sample(xyz, mVal);
677}
678
679template<typename TreeOrAccT>
680__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 2, true>::operator()(const CoordT &ijk) const
681{
682 return ijk == mPos ? mVal[1][1][1] : BaseT::mAcc.getValue(ijk);
683}
684
685template<typename TreeOrAccT>
686template<typename RealT, template<typename...> class Vec3T>
687__hostdev__ bool SampleFromVoxels<TreeOrAccT, 2, true>::zeroCrossing(Vec3T<RealT> xyz) const
688{
689 this->cache(xyz);
690 return BaseT::zeroCrossing(mVal);
691}
692
693template<typename TreeOrAccT>
694template<typename RealT, template<typename...> class Vec3T>
695__hostdev__ void SampleFromVoxels<TreeOrAccT, 2, true>::cache(Vec3T<RealT>& xyz) const
696{
697 CoordT ijk = Floor<CoordT>(xyz);
698 if (ijk != mPos) {
699 mPos = ijk;
700 BaseT::stencil(ijk, mVal);
701 }
702}
703
704template<typename TreeOrAccT>
705template<typename RealT, template<typename...> class Vec3T>
706__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 2, false>::operator()(Vec3T<RealT> xyz) const
707{
708 ValueT val[3][3][3];
709 CoordT ijk = Floor<CoordT>(xyz);
710 BaseT::stencil(ijk, val);
711 return BaseT::sample(xyz, val);
712}
713
714template<typename TreeOrAccT>
715template<typename RealT, template<typename...> class Vec3T>
716__hostdev__ bool SampleFromVoxels<TreeOrAccT, 2, false>::zeroCrossing(Vec3T<RealT> xyz) const
717{
718 ValueT val[3][3][3];
719 CoordT ijk = Floor<CoordT>(xyz);
720 BaseT::stencil(ijk, val);
721 return BaseT::zeroCrossing(val);
722}
723
724// ------------------------------> TricubicSampler <--------------------------------------
725
726/// @brief Tri-cubic sampler, i.e. third order, interpolator.
727///
728/// @details See the following paper for implementation details:
729/// Lekien, F. and Marsden, J.: Tricubic interpolation in three dimensions.
730/// In: International Journal for Numerical Methods
731/// in Engineering (2005), No. 63, p. 455-471
732
733template<typename TreeOrAccT>
734class TricubicSampler
735{
736protected:
737 using ValueT = typename TreeOrAccT::ValueType;
738 using CoordT = typename TreeOrAccT::CoordType;
739
740 const TreeOrAccT& mAcc;
741
742public:
743 /// @brief Construction from a Tree or ReadAccessor
744 __hostdev__ TricubicSampler(const TreeOrAccT& acc)
745 : mAcc(acc)
746 {
747 }
748
749 __hostdev__ const TreeOrAccT& accessor() const { return mAcc; }
750
751 /// @brief Extract the stencil of 8 values
752 inline __hostdev__ void stencil(const CoordT& ijk, ValueT (&c)[64]) const;
753
754 template<typename RealT, template<typename...> class Vec3T>
755 static inline __hostdev__ ValueT sample(const Vec3T<RealT> &uvw, const ValueT (&c)[64]);
756}; // TricubicSampler
757
758template<typename TreeOrAccT>
759__hostdev__ void TricubicSampler<TreeOrAccT>::stencil(const CoordT& ijk, ValueT (&C)[64]) const
760{
761 auto fetch = [&](int i, int j, int k) -> ValueT& { return C[((i + 1) << 4) + ((j + 1) << 2) + k + 1]; };
762
763 // fetch 64 point stencil values
764 for (int i = -1; i < 3; ++i) {
765 for (int j = -1; j < 3; ++j) {
766 fetch(i, j, -1) = mAcc.getValue(ijk + CoordT(i, j, -1));
767 fetch(i, j, 0) = mAcc.getValue(ijk + CoordT(i, j, 0));
768 fetch(i, j, 1) = mAcc.getValue(ijk + CoordT(i, j, 1));
769 fetch(i, j, 2) = mAcc.getValue(ijk + CoordT(i, j, 2));
770 }
771 }
772 // NOTE: these weights truncate to zero for integral ValueT, see detail::toReal
773 const ValueT half = ValueT(0.5), quarter = ValueT(0.25), eighth = ValueT(0.125);
774 const ValueT X[64] = {// values of f(x,y,z) at the 8 corners (each from 1 stencil value).
775 fetch(0, 0, 0),
776 fetch(1, 0, 0),
777 fetch(0, 1, 0),
778 fetch(1, 1, 0),
779 fetch(0, 0, 1),
780 fetch(1, 0, 1),
781 fetch(0, 1, 1),
782 fetch(1, 1, 1),
783 // values of df/dx at the 8 corners (each from 2 stencil values).
784 half * (fetch(1, 0, 0) - fetch(-1, 0, 0)),
785 half * (fetch(2, 0, 0) - fetch(0, 0, 0)),
786 half * (fetch(1, 1, 0) - fetch(-1, 1, 0)),
787 half * (fetch(2, 1, 0) - fetch(0, 1, 0)),
788 half * (fetch(1, 0, 1) - fetch(-1, 0, 1)),
789 half * (fetch(2, 0, 1) - fetch(0, 0, 1)),
790 half * (fetch(1, 1, 1) - fetch(-1, 1, 1)),
791 half * (fetch(2, 1, 1) - fetch(0, 1, 1)),
792 // values of df/dy at the 8 corners (each from 2 stencil values).
793 half * (fetch(0, 1, 0) - fetch(0, -1, 0)),
794 half * (fetch(1, 1, 0) - fetch(1, -1, 0)),
795 half * (fetch(0, 2, 0) - fetch(0, 0, 0)),
796 half * (fetch(1, 2, 0) - fetch(1, 0, 0)),
797 half * (fetch(0, 1, 1) - fetch(0, -1, 1)),
798 half * (fetch(1, 1, 1) - fetch(1, -1, 1)),
799 half * (fetch(0, 2, 1) - fetch(0, 0, 1)),
800 half * (fetch(1, 2, 1) - fetch(1, 0, 1)),
801 // values of df/dz at the 8 corners (each from 2 stencil values).
802 half * (fetch(0, 0, 1) - fetch(0, 0, -1)),
803 half * (fetch(1, 0, 1) - fetch(1, 0, -1)),
804 half * (fetch(0, 1, 1) - fetch(0, 1, -1)),
805 half * (fetch(1, 1, 1) - fetch(1, 1, -1)),
806 half * (fetch(0, 0, 2) - fetch(0, 0, 0)),
807 half * (fetch(1, 0, 2) - fetch(1, 0, 0)),
808 half * (fetch(0, 1, 2) - fetch(0, 1, 0)),
809 half * (fetch(1, 1, 2) - fetch(1, 1, 0)),
810 // values of d2f/dxdy at the 8 corners (each from 4 stencil values).
811 quarter * (fetch(1, 1, 0) - fetch(-1, 1, 0) - fetch(1, -1, 0) + fetch(-1, -1, 0)),
812 quarter * (fetch(2, 1, 0) - fetch(0, 1, 0) - fetch(2, -1, 0) + fetch(0, -1, 0)),
813 quarter * (fetch(1, 2, 0) - fetch(-1, 2, 0) - fetch(1, 0, 0) + fetch(-1, 0, 0)),
814 quarter * (fetch(2, 2, 0) - fetch(0, 2, 0) - fetch(2, 0, 0) + fetch(0, 0, 0)),
815 quarter * (fetch(1, 1, 1) - fetch(-1, 1, 1) - fetch(1, -1, 1) + fetch(-1, -1, 1)),
816 quarter * (fetch(2, 1, 1) - fetch(0, 1, 1) - fetch(2, -1, 1) + fetch(0, -1, 1)),
817 quarter * (fetch(1, 2, 1) - fetch(-1, 2, 1) - fetch(1, 0, 1) + fetch(-1, 0, 1)),
818 quarter * (fetch(2, 2, 1) - fetch(0, 2, 1) - fetch(2, 0, 1) + fetch(0, 0, 1)),
819 // values of d2f/dxdz at the 8 corners (each from 4 stencil values).
820 quarter * (fetch(1, 0, 1) - fetch(-1, 0, 1) - fetch(1, 0, -1) + fetch(-1, 0, -1)),
821 quarter * (fetch(2, 0, 1) - fetch(0, 0, 1) - fetch(2, 0, -1) + fetch(0, 0, -1)),
822 quarter * (fetch(1, 1, 1) - fetch(-1, 1, 1) - fetch(1, 1, -1) + fetch(-1, 1, -1)),
823 quarter * (fetch(2, 1, 1) - fetch(0, 1, 1) - fetch(2, 1, -1) + fetch(0, 1, -1)),
824 quarter * (fetch(1, 0, 2) - fetch(-1, 0, 2) - fetch(1, 0, 0) + fetch(-1, 0, 0)),
825 quarter * (fetch(2, 0, 2) - fetch(0, 0, 2) - fetch(2, 0, 0) + fetch(0, 0, 0)),
826 quarter * (fetch(1, 1, 2) - fetch(-1, 1, 2) - fetch(1, 1, 0) + fetch(-1, 1, 0)),
827 quarter * (fetch(2, 1, 2) - fetch(0, 1, 2) - fetch(2, 1, 0) + fetch(0, 1, 0)),
828 // values of d2f/dydz at the 8 corners (each from 4 stencil values).
829 quarter * (fetch(0, 1, 1) - fetch(0, -1, 1) - fetch(0, 1, -1) + fetch(0, -1, -1)),
830 quarter * (fetch(1, 1, 1) - fetch(1, -1, 1) - fetch(1, 1, -1) + fetch(1, -1, -1)),
831 quarter * (fetch(0, 2, 1) - fetch(0, 0, 1) - fetch(0, 2, -1) + fetch(0, 0, -1)),
832 quarter * (fetch(1, 2, 1) - fetch(1, 0, 1) - fetch(1, 2, -1) + fetch(1, 0, -1)),
833 quarter * (fetch(0, 1, 2) - fetch(0, -1, 2) - fetch(0, 1, 0) + fetch(0, -1, 0)),
834 quarter * (fetch(1, 1, 2) - fetch(1, -1, 2) - fetch(1, 1, 0) + fetch(1, -1, 0)),
835 quarter * (fetch(0, 2, 2) - fetch(0, 0, 2) - fetch(0, 2, 0) + fetch(0, 0, 0)),
836 quarter * (fetch(1, 2, 2) - fetch(1, 0, 2) - fetch(1, 2, 0) + fetch(1, 0, 0)),
837 // values of d3f/dxdydz at the 8 corners (each from 8 stencil values).
838 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)),
839 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)),
840 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)),
841 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)),
842 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)),
843 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)),
844 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)),
845 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))};
846
847 // 4Kb of static table (int8_t has a range of -127 -> 127 which suffices)
848 static const int8_t A[64][64] = {
849 {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},
850 {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},
851 {-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},
852 {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},
853 {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},
854 {0, 0, 0, 0, 0, 0, 0, 0, 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},
855 {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},
856 {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},
857 {-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},
858 {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},
859 {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},
860 {-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},
861 {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},
862 {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},
863 {-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},
864 {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},
865 {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},
866 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
867 {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},
868 {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},
869 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
870 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
871 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
872 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
873 {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},
874 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
875 {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},
876 {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},
877 {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},
878 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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},
879 {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},
880 {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},
881 {-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},
882 {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},
883 {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},
884 {-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},
885 {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},
886 {0, 0, 0, 0, 0, 0, 0, 0, 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},
887 {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},
888 {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},
889 {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},
890 {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},
891 {-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},
892 {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},
893 {-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},
894 {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},
895 {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},
896 {-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},
897 {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},
898 {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},
899 {-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},
900 {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},
901 {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},
902 {0, 0, 0, 0, 0, 0, 0, 0, 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},
903 {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},
904 {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},
905 {-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},
906 {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},
907 {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},
908 {-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},
909 {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},
910 {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},
911 {-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},
912 {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}};
913
914 for (int i = 0; i < 64; ++i) { // C = A * X
915 C[i] = ValueT(0);
916#if 0
917 for (int j = 0; j < 64; j += 4) {
918 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]))));
919 }
920#else
921 for (int j = 0; j < 64; j += 4) {
922 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];
923 }
924#endif
925 }
926}
927
928template<typename TreeOrAccT>
929template<typename RealT, template<typename...> class Vec3T>
930__hostdev__ typename TreeOrAccT::ValueType TricubicSampler<TreeOrAccT>::sample(const Vec3T<RealT> &xyz, const ValueT (&C)[64])
931{
932 ValueT zPow(1), sum(0);
933 for (int k = 0, n = 0; k < 4; ++k) {
934 ValueT yPow(1);
935 for (int j = 0; j < 4; ++j, n += 4) {
936#if 0
937 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);
938#else
939 using detail::toReal;
940 sum = ValueT(toReal<RealT>(sum) + toReal<RealT>(yPow * zPow) * (toReal<RealT>(C[n]) + xyz[0] *
941 (toReal<RealT>(C[n + 1]) + xyz[0] * (toReal<RealT>(C[n + 2]) + xyz[0] * toReal<RealT>(C[n + 3])))));
942#endif
943 yPow = detail::scale(yPow, xyz[1]);
944 }
945 zPow = detail::scale(zPow, xyz[2]);
946 }
947 return sum;
948}
949
950template<typename TreeOrAccT>
951class SampleFromVoxels<TreeOrAccT, 3, true> : public TricubicSampler<TreeOrAccT>
952{
953 using BaseT = TricubicSampler<TreeOrAccT>;
954 using ValueT = typename TreeOrAccT::ValueType;
955 using CoordT = typename TreeOrAccT::CoordType;
956
957 mutable CoordT mPos;
958 mutable ValueT mC[64];
959
960 template<typename RealT, template<typename...> class Vec3T>
961 __hostdev__ void cache(Vec3T<RealT>& xyz) const;
962
963public:
964 /// @brief Construction from a Tree or ReadAccessor
965 __hostdev__ SampleFromVoxels(const TreeOrAccT& acc)
966 : BaseT(acc)
967 {
968 }
969
970 /// @note xyz is in index space space
971 template<typename RealT, template<typename...> class Vec3T>
972 inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
973
974 // @brief Return value at the coordinate @a ijk in index space space
975 __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
976
977}; // SampleFromVoxels<TreeOrAccT, 3, true>
978
979template<typename TreeOrAccT>
980template<typename RealT, template<typename...> class Vec3T>
981__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 3, true>::operator()(Vec3T<RealT> xyz) const
982{
983 this->cache(xyz);
984 return BaseT::sample(xyz, mC);
985}
986
987template<typename TreeOrAccT>
988template<typename RealT, template<typename...> class Vec3T>
989__hostdev__ void SampleFromVoxels<TreeOrAccT, 3, true>::cache(Vec3T<RealT>& xyz) const
990{
991 CoordT ijk = Floor<CoordT>(xyz);
992 if (ijk != mPos) {
993 mPos = ijk;
994 BaseT::stencil(ijk, mC);
995 }
996}
997
998template<typename TreeOrAccT>
999class SampleFromVoxels<TreeOrAccT, 3, false> : public TricubicSampler<TreeOrAccT>
1000{
1001 using BaseT = TricubicSampler<TreeOrAccT>;
1002 using ValueT = typename TreeOrAccT::ValueType;
1003 using CoordT = typename TreeOrAccT::CoordType;
1004
1005public:
1006 /// @brief Construction from a Tree or ReadAccessor
1007 __hostdev__ SampleFromVoxels(const TreeOrAccT& acc)
1008 : BaseT(acc)
1009 {
1010 }
1011
1012 /// @note xyz is in index space space
1013 template<typename RealT, template<typename...> class Vec3T>
1014 inline __hostdev__ ValueT operator()(Vec3T<RealT> xyz) const;
1015
1016 __hostdev__ ValueT operator()(const CoordT &ijk) const {return BaseT::mAcc.getValue(ijk);}
1017
1018}; // SampleFromVoxels<TreeOrAccT, 3, true>
1019
1020template<typename TreeOrAccT>
1021template<typename RealT, template<typename...> class Vec3T>
1022__hostdev__ typename TreeOrAccT::ValueType SampleFromVoxels<TreeOrAccT, 3, false>::operator()(Vec3T<RealT> xyz) const
1023{
1024 ValueT C[64];
1025 CoordT ijk = Floor<CoordT>(xyz);
1026 BaseT::stencil(ijk, C);
1027 return BaseT::sample(xyz, C);
1028}
1029
1030}// namespace math
1031
1032} // namespace nanovdb
1033
1034#endif // NANOVDB_SAMPLE_FROM_VOXELS_H_HAS_BEEN_INCLUDED
__hostdev__
#define __hostdev__
Definition SampleFromVoxels.h:29
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:348
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:341
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:404
nanovdb::math::SampleFromVoxels< TreeOrAccT, 1, true >::SampleFromVoxels
SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:380
nanovdb::math::SampleFromVoxels< TreeOrAccT, 2, false >::operator()
ValueT operator()(const CoordT &ijk) const
Definition SampleFromVoxels.h:624
nanovdb::math::SampleFromVoxels< TreeOrAccT, 2, false >::SampleFromVoxels
SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:618
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:667
nanovdb::math::SampleFromVoxels< TreeOrAccT, 2, true >::SampleFromVoxels
SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:650
nanovdb::math::SampleFromVoxels< TreeOrAccT, 3, false >::operator()
ValueT operator()(const CoordT &ijk) const
Definition SampleFromVoxels.h:1016
nanovdb::math::SampleFromVoxels< TreeOrAccT, 3, false >::SampleFromVoxels
SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:1007
nanovdb::math::SampleFromVoxels< TreeOrAccT, 3, true >::operator()
ValueT operator()(const CoordT &ijk) const
Definition SampleFromVoxels.h:975
nanovdb::math::SampleFromVoxels< TreeOrAccT, 3, true >::SampleFromVoxels
SampleFromVoxels(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:965
nanovdb::math::SampleFromVoxels
Definition SampleFromVoxels.h:40
nanovdb::math::TricubicSampler::accessor
const TreeOrAccT & accessor() const
Definition SampleFromVoxels.h:749
nanovdb::math::TricubicSampler::TricubicSampler
TricubicSampler(const TreeOrAccT &acc)
Construction from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:744
nanovdb::math::TricubicSampler::mAcc
const TreeOrAccT & mAcc
Definition SampleFromVoxels.h:740
nanovdb::math::TricubicSampler::CoordT
typename TreeOrAccT::CoordType CoordT
Definition SampleFromVoxels.h:738
nanovdb::math::TricubicSampler::stencil
void stencil(const CoordT &ijk, ValueT(&c)[64]) const
Extract the stencil of 8 values.
Definition SampleFromVoxels.h:759
nanovdb::math::TricubicSampler::ValueT
typename TreeOrAccT::ValueType ValueT
Definition SampleFromVoxels.h:737
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:241
nanovdb::math::TrilinearSampler::TrilinearSampler
TrilinearSampler(const TreeOrAccT &acc)
Protected constructor from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:224
nanovdb::math::TrilinearSampler::accessor
const TreeOrAccT & accessor() const
Definition SampleFromVoxels.h:226
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:216
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:221
nanovdb::math::TrilinearSampler::zeroCrossing
static bool zeroCrossing(const ValueT(&v)[2][2][2])
Definition SampleFromVoxels.h:317
nanovdb::math::TrilinearSampler::CoordT
typename TreeOrAccT::CoordType CoordT
Definition SampleFromVoxels.h:220
nanovdb::math::TrilinearSampler::ValueT
typename TreeOrAccT::ValueType ValueT
Definition SampleFromVoxels.h:219
nanovdb::math::TriquadraticSampler::TriquadraticSampler
TriquadraticSampler(const TreeOrAccT &acc)
Protected constructor from a Tree or ReadAccessor.
Definition SampleFromVoxels.h:544
nanovdb::math::TriquadraticSampler::stencil
void stencil(const CoordT &ijk, ValueT(&v)[3][3][3]) const
Extract the stencil of 27 values.
Definition SampleFromVoxels.h:558
nanovdb::math::TriquadraticSampler::accessor
const TreeOrAccT & accessor() const
Definition SampleFromVoxels.h:546
nanovdb::math::TriquadraticSampler::mAcc
const TreeOrAccT & mAcc
Definition SampleFromVoxels.h:536
nanovdb::math::TriquadraticSampler::ORDER
static const int ORDER
Definition SampleFromVoxels.h:541
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:540
nanovdb::math::TriquadraticSampler::ValueT
typename TreeOrAccT::ValueType ValueT
Definition SampleFromVoxels.h:539
nanovdb::math::TriquadraticSampler::zeroCrossing
static bool zeroCrossing(const ValueT(&v)[3][3][3])
Definition SampleFromVoxels.h:594
detail
Definition VoxToNanoVDB.h:15
nanovdb::math::detail::toReal
constexpr auto toReal(const T &v)
Explicitly convert an integral v to RealT, i.e. the implicit promotion in mixed integer/floating-poin...
Definition SampleFromVoxels.h:187
nanovdb::math::detail::scale
T scale(const T &v, S s)
Preserve NanoVDB vectors' component-precision scalar multiplication. For other types,...
Definition SampleFromVoxels.h:199
nanovdb::math::Floor
__hostdev__ int32_t Floor(float x) noexcept
Floor of x as a 32-bit signed integer (truncation toward -inf).
Definition Math.h:234
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::Round
__hostdev__ int32_t Round(float x) noexcept
Round x to the nearest 32-bit signed integer, with halfway cases rounded toward +inf.
Definition Math.h:308
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:194
nanovdb
Defines a simple memory pool used to call cub functions that use dynamic temporary storage.
Definition GridHandle.h:34
Math.h
Math functions and classes.
nanovdb::math::detail::MakeVoid
Scalar component type of a scalar or vector value type, e.g. float for both float and Vec3<float>
Definition SampleFromVoxels.h:175
nanovdb::math::detail::MakeVoid::type
void type
Definition SampleFromVoxels.h:175
nanovdb::math::detail::ScalarOf
Definition SampleFromVoxels.h:177
nanovdb::math::detail::ScalarOf::type
T type
Definition SampleFromVoxels.h:177
nanovdb::util::is_floating_point::value
static constexpr bool value
Definition Util.h:375
nanovdb::util::is_specialization::value
static const bool value
Definition Util.h:527