OpenVDB 10.0.1
Loading...
Searching...
No Matches
Interpolation.h
Go to the documentation of this file.
1// Copyright Contributors to the OpenVDB Project
2// SPDX-License-Identifier: MPL-2.0
3
4/// @file Interpolation.h
5///
6/// Sampler classes such as PointSampler and BoxSampler that are intended for use
7/// with tools::GridTransformer should operate in voxel space and must adhere to
8/// the interface described in the example below:
9/// @code
10/// struct MySampler
11/// {
12/// // Return a short name that can be used to identify this sampler
13/// // in error messages and elsewhere.
14/// const char* name() { return "mysampler"; }
15///
16/// // Return the radius of the sampling kernel in voxels, not including
17/// // the center voxel. This is the number of voxels of padding that
18/// // are added to all sides of a volume as a result of resampling.
19/// int radius() { return 2; }
20///
21/// // Return true if scaling by a factor smaller than 0.5 (along any axis)
22/// // should be handled via a mipmapping-like scheme of successive halvings
23/// // of a grid's resolution, until the remaining scale factor is
24/// // greater than or equal to 1/2. Set this to false only when high-quality
25/// // scaling is not required.
26/// bool mipmap() { return true; }
27///
28/// // Specify if sampling at a location that is collocated with a grid point
29/// // is guaranteed to return the exact value at that grid point.
30/// // For most sampling kernels, this should be false.
31/// bool consistent() { return false; }
32///
33/// // Sample the tree at the given coordinates and return the result in val.
34/// // Return true if the sampled value is active.
35/// template<class TreeT>
36/// bool sample(const TreeT& tree, const Vec3R& coord, typename TreeT::ValueType& val);
37/// };
38/// @endcode
39
40#ifndef OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
41#define OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
42
43#include <openvdb/version.h> // for OPENVDB_VERSION_NAME
44#include <openvdb/Platform.h> // for round()
45#include <openvdb/math/Math.h>// for SmoothUnitStep
46#include <openvdb/math/Transform.h> // for Transform
47#include <openvdb/Grid.h>
49#include <cmath>
50#include <type_traits>
51
52namespace openvdb {
54namespace OPENVDB_VERSION_NAME {
55namespace tools {
56
57/// @brief Provises a unified interface for sampling, i.e. interpolation.
58/// @details Order = 0: closest point
59/// Order = 1: tri-linear
60/// Order = 2: tri-quadratic
61/// Staggered: Set to true for MAC grids
62template <size_t Order, bool Staggered = false>
63struct Sampler
64{
65 static_assert(Order < 3, "Samplers of order higher than 2 are not supported");
66 static const char* name();
67 static int radius();
68 static bool mipmap();
69 static bool consistent();
70 static bool staggered();
71 static size_t order();
72
73 /// @brief Sample @a inTree at the floating-point index coordinate @a inCoord
74 /// and store the result in @a result.
75 ///
76 /// @return @c true if the sampled value is active.
77 template<class TreeT>
78 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
79 typename TreeT::ValueType& result);
80
81 /// @brief Sample @a inTree at the floating-point index coordinate @a inCoord.
82 ///
83 /// @return the reconstructed value
84 template<class TreeT>
85 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
86};
87
88//////////////////////////////////////// Non-Staggered Samplers
89
90// The following samplers operate in voxel space.
91// When the samplers are applied to grids holding vector or other non-scalar data,
92// the data is assumed to be collocated. For example, using the BoxSampler on a grid
93// with ValueType Vec3f assumes that all three elements in a vector can be assigned
94// the same physical location. Consider using the GridSampler below instead.
95
97{
98 static const char* name() { return "point"; }
99 static int radius() { return 0; }
100 static bool mipmap() { return false; }
101 static bool consistent() { return true; }
102 static bool staggered() { return false; }
103 static size_t order() { return 0; }
104
105 /// @brief Sample @a inTree at the nearest neighbor to @a inCoord
106 /// and store the result in @a result.
107 /// @return @c true if the sampled value is active.
108 template<class TreeT>
109 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
110 typename TreeT::ValueType& result);
111
112 /// @brief Sample @a inTree at the nearest neighbor to @a inCoord
113 /// @return the reconstructed value
114 template<class TreeT>
115 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
116};
117
118
120{
121 static const char* name() { return "box"; }
122 static int radius() { return 1; }
123 static bool mipmap() { return true; }
124 static bool consistent() { return true; }
125 static bool staggered() { return false; }
126 static size_t order() { return 1; }
127
128 /// @brief Trilinearly reconstruct @a inTree at @a inCoord
129 /// and store the result in @a result.
130 /// @return @c true if any one of the sampled values is active.
131 template<class TreeT>
132 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
133 typename TreeT::ValueType& result);
134
135 /// @brief Trilinearly reconstruct @a inTree at @a inCoord.
136 /// @return the reconstructed value
137 template<class TreeT>
138 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
139
140 /// @brief Import all eight values from @a inTree to support
141 /// tri-linear interpolation.
142 template<class ValueT, class TreeT, size_t N>
143 static inline void getValues(ValueT (&data)[N][N][N], const TreeT& inTree, Coord ijk);
144
145 /// @brief Import all eight values from @a inTree to support
146 /// tri-linear interpolation.
147 /// @return @c true if any of the eight values are active
148 template<class ValueT, class TreeT, size_t N>
149 static inline bool probeValues(ValueT (&data)[N][N][N], const TreeT& inTree, Coord ijk);
150
151 /// @brief Find the minimum and maximum values of the eight cell
152 /// values in @ data.
153 template<class ValueT, size_t N>
154 static inline void extrema(ValueT (&data)[N][N][N], ValueT& vMin, ValueT& vMax);
155
156 /// @return the tri-linear interpolation with the unit cell coordinates @a uvw
157 template<class ValueT, size_t N>
158 static inline ValueT trilinearInterpolation(ValueT (&data)[N][N][N], const Vec3R& uvw);
159};
160
161
163{
164 static const char* name() { return "quadratic"; }
165 static int radius() { return 1; }
166 static bool mipmap() { return true; }
167 static bool consistent() { return false; }
168 static bool staggered() { return false; }
169 static size_t order() { return 2; }
170
171 /// @brief Triquadratically reconstruct @a inTree at @a inCoord
172 /// and store the result in @a result.
173 /// @return @c true if any one of the sampled values is active.
174 template<class TreeT>
175 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
176 typename TreeT::ValueType& result);
177
178 /// @brief Triquadratically reconstruct @a inTree at to @a inCoord.
179 /// @return the reconstructed value
180 template<class TreeT>
181 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
182
183 template<class ValueT, size_t N>
184 static inline ValueT triquadraticInterpolation(ValueT (&data)[N][N][N], const Vec3R& uvw);
185};
186
187
188//////////////////////////////////////// Staggered Samplers
189
190
191// The following samplers operate in voxel space and are designed for Vec3
192// staggered grid data (e.g., fluid simulations using the Marker-and-Cell approach
193// associate elements of the velocity vector with different physical locations:
194// the faces of a cube).
195
197{
198 static const char* name() { return "point"; }
199 static int radius() { return 0; }
200 static bool mipmap() { return false; }
201 static bool consistent() { return false; }
202 static bool staggered() { return true; }
203 static size_t order() { return 0; }
204
205 /// @brief Sample @a inTree at the nearest neighbor to @a inCoord
206 /// and store the result in @a result.
207 /// @return true if the sampled value is active.
208 template<class TreeT>
209 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
210 typename TreeT::ValueType& result);
211
212 /// @brief Sample @a inTree at the nearest neighbor to @a inCoord
213 /// @return the reconstructed value
214 template<class TreeT>
215 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
216};
217
218
220{
221 static const char* name() { return "box"; }
222 static int radius() { return 1; }
223 static bool mipmap() { return true; }
224 static bool consistent() { return false; }
225 static bool staggered() { return true; }
226 static size_t order() { return 1; }
227
228 /// @brief Trilinearly reconstruct @a inTree at @a inCoord
229 /// and store the result in @a result.
230 /// @return true if any one of the sampled value is active.
231 template<class TreeT>
232 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
233 typename TreeT::ValueType& result);
234
235 /// @brief Trilinearly reconstruct @a inTree at @a inCoord.
236 /// @return the reconstructed value
237 template<class TreeT>
238 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
239};
240
241
243{
244 static const char* name() { return "quadratic"; }
245 static int radius() { return 1; }
246 static bool mipmap() { return true; }
247 static bool consistent() { return false; }
248 static bool staggered() { return true; }
249 static size_t order() { return 2; }
250
251 /// @brief Triquadratically reconstruct @a inTree at @a inCoord
252 /// and store the result in @a result.
253 /// @return true if any one of the sampled values is active.
254 template<class TreeT>
255 static bool sample(const TreeT& inTree, const Vec3R& inCoord,
256 typename TreeT::ValueType& result);
257
258 /// @brief Triquadratically reconstruct @a inTree at to @a inCoord.
259 /// @return the reconstructed value
260 template<class TreeT>
261 static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
262};
263
264
265//////////////////////////////////////// GridSampler
266
267
268/// @brief Class that provides the interface for continuous sampling
269/// of values in a tree.
270///
271/// @details Since trees support only discrete voxel sampling, TreeSampler
272/// must be used to sample arbitrary continuous points in (world or
273/// index) space.
274///
275/// @warning This implementation of the GridSampler stores a pointer
276/// to a Tree for value access. While this is thread-safe it is
277/// uncached and hence slow compared to using a
278/// ValueAccessor. Consequently it is normally advisable to use the
279/// template specialization below that employs a
280/// ValueAccessor. However, care must be taken when dealing with
281/// multi-threading (see warning below).
282template<typename GridOrTreeType, typename SamplerType>
284{
285public:
287 using ValueType = typename GridOrTreeType::ValueType;
291
292 /// @param grid a grid to be sampled
293 explicit GridSampler(const GridType& grid)
294 : mTree(&(grid.tree())), mTransform(&(grid.transform())) {}
295
296 /// @param tree a tree to be sampled, or a ValueAccessor for the tree
297 /// @param transform is used when sampling world space locations.
298 GridSampler(const TreeType& tree, const math::Transform& transform)
299 : mTree(&tree), mTransform(&transform) {}
300
301 const math::Transform& transform() const { return *mTransform; }
302
303 /// @brief Sample a point in index space in the grid.
304 /// @param x Fractional x-coordinate of point in index-coordinates of grid
305 /// @param y Fractional y-coordinate of point in index-coordinates of grid
306 /// @param z Fractional z-coordinate of point in index-coordinates of grid
307 template<typename RealType>
308 ValueType sampleVoxel(const RealType& x, const RealType& y, const RealType& z) const
309 {
310 return this->isSample(Vec3d(x,y,z));
311 }
312
313 /// @brief Sample value in integer index space
314 /// @param i Integer x-coordinate in index space
315 /// @param j Integer y-coordinate in index space
316 /// @param k Integer x-coordinate in index space
318 typename Coord::ValueType j,
319 typename Coord::ValueType k) const
320 {
321 return this->isSample(Coord(i,j,k));
322 }
323
324 /// @brief Sample value in integer index space
325 /// @param ijk the location in index space
326 ValueType isSample(const Coord& ijk) const { return mTree->getValue(ijk); }
327
328 /// @brief Sample in fractional index space
329 /// @param ispoint the location in index space
330 ValueType isSample(const Vec3d& ispoint) const
331 {
332 ValueType result = zeroVal<ValueType>();
333 SamplerType::sample(*mTree, ispoint, result);
334 return result;
335 }
336
337 /// @brief Sample in world space
338 /// @param wspoint the location in world space
339 ValueType wsSample(const Vec3d& wspoint) const
340 {
341 ValueType result = zeroVal<ValueType>();
342 SamplerType::sample(*mTree, mTransform->worldToIndex(wspoint), result);
343 return result;
344 }
345
346private:
347 const TreeType* mTree;
348 const math::Transform* mTransform;
349}; // class GridSampler
350
351
352/// @brief Specialization of GridSampler for construction from a ValueAccessor type
353///
354/// @note This version should normally be favored over the one above
355/// that takes a Grid or Tree. The reason is this version uses a
356/// ValueAccessor that performs fast (cached) access where the
357/// tree-based flavor performs slower (uncached) access.
358///
359/// @warning Since this version stores a pointer to an (externally
360/// allocated) value accessor it is not threadsafe. Hence each thread
361/// should have its own instance of a GridSampler constructed from a
362/// local ValueAccessor. Alternatively the Grid/Tree-based GridSampler
363/// is threadsafe, but also slower.
364template<typename TreeT, typename SamplerType>
365class GridSampler<tree::ValueAccessor<TreeT>, SamplerType>
366{
367public:
369 using ValueType = typename TreeT::ValueType;
370 using TreeType = TreeT;
373
374 /// @param acc a ValueAccessor to be sampled
375 /// @param transform is used when sampling world space locations.
377 const math::Transform& transform)
378 : mAccessor(&acc), mTransform(&transform) {}
379
380 const math::Transform& transform() const { return *mTransform; }
381
382 /// @brief Sample a point in index space in the grid.
383 /// @param x Fractional x-coordinate of point in index-coordinates of grid
384 /// @param y Fractional y-coordinate of point in index-coordinates of grid
385 /// @param z Fractional z-coordinate of point in index-coordinates of grid
386 template<typename RealType>
387 ValueType sampleVoxel(const RealType& x, const RealType& y, const RealType& z) const
388 {
389 return this->isSample(Vec3d(x,y,z));
390 }
391
392 /// @brief Sample value in integer index space
393 /// @param i Integer x-coordinate in index space
394 /// @param j Integer y-coordinate in index space
395 /// @param k Integer x-coordinate in index space
397 typename Coord::ValueType j,
398 typename Coord::ValueType k) const
399 {
400 return this->isSample(Coord(i,j,k));
401 }
402
403 /// @brief Sample value in integer index space
404 /// @param ijk the location in index space
405 ValueType isSample(const Coord& ijk) const { return mAccessor->getValue(ijk); }
406
407 /// @brief Sample in fractional index space
408 /// @param ispoint the location in index space
409 ValueType isSample(const Vec3d& ispoint) const
410 {
411 ValueType result = zeroVal<ValueType>();
412 SamplerType::sample(*mAccessor, ispoint, result);
413 return result;
414 }
415
416 /// @brief Sample in world space
417 /// @param wspoint the location in world space
418 ValueType wsSample(const Vec3d& wspoint) const
419 {
420 ValueType result = zeroVal<ValueType>();
421 SamplerType::sample(*mAccessor, mTransform->worldToIndex(wspoint), result);
422 return result;
423 }
424
425private:
426 const AccessorType* mAccessor;//not thread-safe!
427 const math::Transform* mTransform;
428};//Specialization of GridSampler
429
430
431//////////////////////////////////////// DualGridSampler
432
433
434/// @brief This is a simple convenience class that allows for sampling
435/// from a source grid into the index space of a target grid. At
436/// construction the source and target grids are checked for alignment
437/// which potentially renders interpolation unnecessary. Else
438/// interpolation is performed according to the templated Sampler
439/// type.
440///
441/// @warning For performance reasons the check for alignment of the
442/// two grids is only performed at construction time!
443template<typename GridOrTreeT,
444 typename SamplerT>
446{
447public:
448 using ValueType = typename GridOrTreeT::ValueType;
452
453 /// @brief Grid and transform constructor.
454 /// @param sourceGrid Source grid.
455 /// @param targetXform Transform of the target grid.
456 DualGridSampler(const GridType& sourceGrid,
457 const math::Transform& targetXform)
458 : mSourceTree(&(sourceGrid.tree()))
459 , mSourceXform(&(sourceGrid.transform()))
460 , mTargetXform(&targetXform)
461 , mAligned(targetXform == *mSourceXform)
462 {
463 }
464 /// @brief Tree and transform constructor.
465 /// @param sourceTree Source tree.
466 /// @param sourceXform Transform of the source grid.
467 /// @param targetXform Transform of the target grid.
468 DualGridSampler(const TreeType& sourceTree,
469 const math::Transform& sourceXform,
470 const math::Transform& targetXform)
471 : mSourceTree(&sourceTree)
472 , mSourceXform(&sourceXform)
473 , mTargetXform(&targetXform)
474 , mAligned(targetXform == sourceXform)
475 {
476 }
477 /// @brief Return the value of the source grid at the index
478 /// coordinates, ijk, relative to the target grid (or its tranform).
479 inline ValueType operator()(const Coord& ijk) const
480 {
481 if (mAligned) return mSourceTree->getValue(ijk);
482 const Vec3R world = mTargetXform->indexToWorld(ijk);
483 return SamplerT::sample(*mSourceTree, mSourceXform->worldToIndex(world));
484 }
485 /// @brief Return true if the two grids are aligned.
486 inline bool isAligned() const { return mAligned; }
487private:
488 const TreeType* mSourceTree;
489 const math::Transform* mSourceXform;
490 const math::Transform* mTargetXform;
491 const bool mAligned;
492};// DualGridSampler
493
494/// @brief Specialization of DualGridSampler for construction from a ValueAccessor type.
495template<typename TreeT,
496 typename SamplerT>
497class DualGridSampler<tree::ValueAccessor<TreeT>, SamplerT>
498{
499 public:
500 using ValueType = typename TreeT::ValueType;
501 using TreeType = TreeT;
504
505 /// @brief ValueAccessor and transform constructor.
506 /// @param sourceAccessor ValueAccessor into the source grid.
507 /// @param sourceXform Transform for the source grid.
508 /// @param targetXform Transform for the target grid.
509 DualGridSampler(const AccessorType& sourceAccessor,
510 const math::Transform& sourceXform,
511 const math::Transform& targetXform)
512 : mSourceAcc(&sourceAccessor)
513 , mSourceXform(&sourceXform)
514 , mTargetXform(&targetXform)
515 , mAligned(targetXform == sourceXform)
516 {
517 }
518 /// @brief Return the value of the source grid at the index
519 /// coordinates, ijk, relative to the target grid.
520 inline ValueType operator()(const Coord& ijk) const
521 {
522 if (mAligned) return mSourceAcc->getValue(ijk);
523 const Vec3R world = mTargetXform->indexToWorld(ijk);
524 return SamplerT::sample(*mSourceAcc, mSourceXform->worldToIndex(world));
525 }
526 /// @brief Return true if the two grids are aligned.
527 inline bool isAligned() const { return mAligned; }
528private:
529 const AccessorType* mSourceAcc;
530 const math::Transform* mSourceXform;
531 const math::Transform* mTargetXform;
532 const bool mAligned;
533};//Specialization of DualGridSampler
534
535//////////////////////////////////////// AlphaMask
536
537
538// Class to derive the normalized alpha mask
539template <typename GridT,
540 typename MaskT,
541 typename SamplerT = tools::BoxSampler,
542 typename FloatT = float>
544{
545public:
546 static_assert(std::is_floating_point<FloatT>::value,
547 "AlphaMask requires a floating-point value type");
548 using GridType = GridT;
549 using MaskType = MaskT;
550 using SamlerType = SamplerT;
551 using FloatType = FloatT;
552
553 AlphaMask(const GridT& grid, const MaskT& mask, FloatT min, FloatT max, bool invert)
554 : mAcc(mask.tree())
555 , mSampler(mAcc, mask.transform() , grid.transform())
556 , mMin(min)
557 , mInvNorm(1/(max-min))
558 , mInvert(invert)
559 {
560 assert(min < max);
561 }
562
563 inline bool operator()(const Coord& xyz, FloatT& a, FloatT& b) const
564 {
565 a = math::SmoothUnitStep( (mSampler(xyz) - mMin) * mInvNorm );//smooth mapping to 0->1
566 b = 1 - a;
567 if (mInvert) std::swap(a,b);
568 return a>0;
569 }
570
571protected:
572 using AccT = typename MaskType::ConstAccessor;
575 const FloatT mMin, mInvNorm;
576 const bool mInvert;
577};// AlphaMask
578
579////////////////////////////////////////
580
581namespace local_util {
582
583inline Vec3i
585{
586 return Vec3i(int(std::floor(v(0))), int(std::floor(v(1))), int(std::floor(v(2))));
587}
588
589
590inline Vec3i
592{
593 return Vec3i(int(std::ceil(v(0))), int(std::ceil(v(1))), int(std::ceil(v(2))));
594}
595
596
597inline Vec3i
599{
600 return Vec3i(int(::round(v(0))), int(::round(v(1))), int(::round(v(2))));
601}
602
603} // namespace local_util
604
605
606//////////////////////////////////////// PointSampler
607
608
609template<class TreeT>
610inline bool
611PointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
612 typename TreeT::ValueType& result)
613{
614 return inTree.probeValue(Coord(local_util::roundVec3(inCoord)), result);
615}
616
617template<class TreeT>
618inline typename TreeT::ValueType
619PointSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
620{
621 return inTree.getValue(Coord(local_util::roundVec3(inCoord)));
622}
623
624
625//////////////////////////////////////// BoxSampler
626
627template<class ValueT, class TreeT, size_t N>
628inline void
629BoxSampler::getValues(ValueT (&data)[N][N][N], const TreeT& inTree, Coord ijk)
630{
631 data[0][0][0] = inTree.getValue(ijk); // i, j, k
632
633 ijk[2] += 1;
634 data[0][0][1] = inTree.getValue(ijk); // i, j, k + 1
635
636 ijk[1] += 1;
637 data[0][1][1] = inTree.getValue(ijk); // i, j+1, k + 1
638
639 ijk[2] -= 1;
640 data[0][1][0] = inTree.getValue(ijk); // i, j+1, k
641
642 ijk[0] += 1;
643 ijk[1] -= 1;
644 data[1][0][0] = inTree.getValue(ijk); // i+1, j, k
645
646 ijk[2] += 1;
647 data[1][0][1] = inTree.getValue(ijk); // i+1, j, k + 1
648
649 ijk[1] += 1;
650 data[1][1][1] = inTree.getValue(ijk); // i+1, j+1, k + 1
651
652 ijk[2] -= 1;
653 data[1][1][0] = inTree.getValue(ijk); // i+1, j+1, k
654}
655
656template<class ValueT, class TreeT, size_t N>
657inline bool
658BoxSampler::probeValues(ValueT (&data)[N][N][N], const TreeT& inTree, Coord ijk)
659{
660 bool hasActiveValues = false;
661 hasActiveValues |= inTree.probeValue(ijk, data[0][0][0]); // i, j, k
662
663 ijk[2] += 1;
664 hasActiveValues |= inTree.probeValue(ijk, data[0][0][1]); // i, j, k + 1
665
666 ijk[1] += 1;
667 hasActiveValues |= inTree.probeValue(ijk, data[0][1][1]); // i, j+1, k + 1
668
669 ijk[2] -= 1;
670 hasActiveValues |= inTree.probeValue(ijk, data[0][1][0]); // i, j+1, k
671
672 ijk[0] += 1;
673 ijk[1] -= 1;
674 hasActiveValues |= inTree.probeValue(ijk, data[1][0][0]); // i+1, j, k
675
676 ijk[2] += 1;
677 hasActiveValues |= inTree.probeValue(ijk, data[1][0][1]); // i+1, j, k + 1
678
679 ijk[1] += 1;
680 hasActiveValues |= inTree.probeValue(ijk, data[1][1][1]); // i+1, j+1, k + 1
681
682 ijk[2] -= 1;
683 hasActiveValues |= inTree.probeValue(ijk, data[1][1][0]); // i+1, j+1, k
684
685 return hasActiveValues;
686}
687
688template<class ValueT, size_t N>
689inline void
690BoxSampler::extrema(ValueT (&data)[N][N][N], ValueT& vMin, ValueT &vMax)
691{
692 vMin = vMax = data[0][0][0];
693 vMin = math::Min(vMin, data[0][0][1]);
694 vMax = math::Max(vMax, data[0][0][1]);
695 vMin = math::Min(vMin, data[0][1][0]);
696 vMax = math::Max(vMax, data[0][1][0]);
697 vMin = math::Min(vMin, data[0][1][1]);
698 vMax = math::Max(vMax, data[0][1][1]);
699 vMin = math::Min(vMin, data[1][0][0]);
700 vMax = math::Max(vMax, data[1][0][0]);
701 vMin = math::Min(vMin, data[1][0][1]);
702 vMax = math::Max(vMax, data[1][0][1]);
703 vMin = math::Min(vMin, data[1][1][0]);
704 vMax = math::Max(vMax, data[1][1][0]);
705 vMin = math::Min(vMin, data[1][1][1]);
706 vMax = math::Max(vMax, data[1][1][1]);
707}
708
709
710template<class ValueT, size_t N>
711inline ValueT
712BoxSampler::trilinearInterpolation(ValueT (&data)[N][N][N], const Vec3R& uvw)
713{
714 auto _interpolate = [](const ValueT& a, const ValueT& b, double weight)
715 {
717 const auto temp = (b - a) * weight;
719 return static_cast<ValueT>(a + ValueT(temp));
720 };
721
722 // Trilinear interpolation:
723 // The eight surrounding latice values are used to construct the result. \n
724 // result(x,y,z) =
725 // v000 (1-x)(1-y)(1-z) + v001 (1-x)(1-y)z + v010 (1-x)y(1-z) + v011 (1-x)yz
726 // + v100 x(1-y)(1-z) + v101 x(1-y)z + v110 xy(1-z) + v111 xyz
727
728 return _interpolate(
729 _interpolate(
730 _interpolate(data[0][0][0], data[0][0][1], uvw[2]),
731 _interpolate(data[0][1][0], data[0][1][1], uvw[2]),
732 uvw[1]),
733 _interpolate(
734 _interpolate(data[1][0][0], data[1][0][1], uvw[2]),
735 _interpolate(data[1][1][0], data[1][1][1], uvw[2]),
736 uvw[1]),
737 uvw[0]);
738}
739
740
741template<class TreeT>
742inline bool
743BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
744 typename TreeT::ValueType& result)
745{
746 using ValueT = typename TreeT::ValueType;
747
748 const Vec3i inIdx = local_util::floorVec3(inCoord);
749 const Vec3R uvw = inCoord - inIdx;
750
751 // Retrieve the values of the eight voxels surrounding the
752 // fractional source coordinates.
753 ValueT data[2][2][2];
754
755 const bool hasActiveValues = BoxSampler::probeValues(data, inTree, Coord(inIdx));
756
757 result = BoxSampler::trilinearInterpolation(data, uvw);
758
759 return hasActiveValues;
760}
761
762
763template<class TreeT>
764inline typename TreeT::ValueType
765BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
766{
767 using ValueT = typename TreeT::ValueType;
768
769 const Vec3i inIdx = local_util::floorVec3(inCoord);
770 const Vec3R uvw = inCoord - inIdx;
771
772 // Retrieve the values of the eight voxels surrounding the
773 // fractional source coordinates.
774 ValueT data[2][2][2];
775
776 BoxSampler::getValues(data, inTree, Coord(inIdx));
777
778 return BoxSampler::trilinearInterpolation(data, uvw);
779}
780
781
782//////////////////////////////////////// QuadraticSampler
783
784template<class ValueT, size_t N>
785inline ValueT
786QuadraticSampler::triquadraticInterpolation(ValueT (&data)[N][N][N], const Vec3R& uvw)
787{
788 auto _interpolate = [](const ValueT* value, double weight)
789 {
791 const ValueT
792 a = static_cast<ValueT>(0.5 * (value[0] + value[2]) - value[1]),
793 b = static_cast<ValueT>(0.5 * (value[2] - value[0])),
794 c = static_cast<ValueT>(value[1]);
795 const auto temp = weight * (weight * a + b) + c;
797 return static_cast<ValueT>(temp);
798 };
799
800 /// @todo For vector types, interpolate over each component independently.
801 ValueT vx[3];
802 for (int dx = 0; dx < 3; ++dx) {
803 ValueT vy[3];
804 for (int dy = 0; dy < 3; ++dy) {
805 // Fit a parabola to three contiguous samples in z
806 // (at z=-1, z=0 and z=1), then evaluate the parabola at z',
807 // where z' is the fractional part of inCoord.z, i.e.,
808 // inCoord.z - inIdx.z. The coefficients come from solving
809 //
810 // | (-1)^2 -1 1 || a | | v0 |
811 // | 0 0 1 || b | = | v1 |
812 // | 1^2 1 1 || c | | v2 |
813 //
814 // for a, b and c.
815 const ValueT* vz = &data[dx][dy][0];
816 vy[dy] = _interpolate(vz, uvw.z());
817 }//loop over y
818 // Fit a parabola to three interpolated samples in y, then
819 // evaluate the parabola at y', where y' is the fractional
820 // part of inCoord.y.
821 vx[dx] = _interpolate(vy, uvw.y());
822 }//loop over x
823 // Fit a parabola to three interpolated samples in x, then
824 // evaluate the parabola at the fractional part of inCoord.x.
825 return _interpolate(vx, uvw.x());
826}
827
828template<class TreeT>
829inline bool
830QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
831 typename TreeT::ValueType& result)
832{
833 using ValueT = typename TreeT::ValueType;
834
835 const Vec3i inIdx = local_util::floorVec3(inCoord), inLoIdx = inIdx - Vec3i(1, 1, 1);
836 const Vec3R uvw = inCoord - inIdx;
837
838 // Retrieve the values of the 27 voxels surrounding the
839 // fractional source coordinates.
840 bool active = false;
841 ValueT data[3][3][3];
842 for (int dx = 0, ix = inLoIdx.x(); dx < 3; ++dx, ++ix) {
843 for (int dy = 0, iy = inLoIdx.y(); dy < 3; ++dy, ++iy) {
844 for (int dz = 0, iz = inLoIdx.z(); dz < 3; ++dz, ++iz) {
845 if (inTree.probeValue(Coord(ix, iy, iz), data[dx][dy][dz])) active = true;
846 }
847 }
848 }
849
850 result = QuadraticSampler::triquadraticInterpolation(data, uvw);
851
852 return active;
853}
854
855template<class TreeT>
856inline typename TreeT::ValueType
857QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
858{
859 using ValueT = typename TreeT::ValueType;
860
861 const Vec3i inIdx = local_util::floorVec3(inCoord), inLoIdx = inIdx - Vec3i(1, 1, 1);
862 const Vec3R uvw = inCoord - inIdx;
863
864 // Retrieve the values of the 27 voxels surrounding the
865 // fractional source coordinates.
866 ValueT data[3][3][3];
867 for (int dx = 0, ix = inLoIdx.x(); dx < 3; ++dx, ++ix) {
868 for (int dy = 0, iy = inLoIdx.y(); dy < 3; ++dy, ++iy) {
869 for (int dz = 0, iz = inLoIdx.z(); dz < 3; ++dz, ++iz) {
870 data[dx][dy][dz] = inTree.getValue(Coord(ix, iy, iz));
871 }
872 }
873 }
874
875 return QuadraticSampler::triquadraticInterpolation(data, uvw);
876}
877
878
879//////////////////////////////////////// StaggeredPointSampler
880
881
882template<class TreeT>
883inline bool
884StaggeredPointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
885 typename TreeT::ValueType& result)
886{
887 using ValueType = typename TreeT::ValueType;
888
889 ValueType tempX, tempY, tempZ;
890 bool active = false;
891
892 active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
893 active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
894 active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
895
896 result.x() = tempX.x();
897 result.y() = tempY.y();
898 result.z() = tempZ.z();
899
900 return active;
901}
902
903template<class TreeT>
904inline typename TreeT::ValueType
905StaggeredPointSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
906{
907 using ValueT = typename TreeT::ValueType;
908
909 const ValueT tempX = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0.0, 0.0));
910 const ValueT tempY = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.0, 0.5, 0.0));
911 const ValueT tempZ = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.0, 0.0, 0.5));
912
913 return ValueT(tempX.x(), tempY.y(), tempZ.z());
914}
915
916
917//////////////////////////////////////// StaggeredBoxSampler
918
919
920template<class TreeT>
921inline bool
922StaggeredBoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
923 typename TreeT::ValueType& result)
924{
925 using ValueType = typename TreeT::ValueType;
926
927 ValueType tempX, tempY, tempZ;
928 tempX = tempY = tempZ = zeroVal<ValueType>();
929 bool active = false;
930
931 active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
932 active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
933 active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
934
935 result.x() = tempX.x();
936 result.y() = tempY.y();
937 result.z() = tempZ.z();
938
939 return active;
940}
941
942template<class TreeT>
943inline typename TreeT::ValueType
944StaggeredBoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
945{
946 using ValueT = typename TreeT::ValueType;
947
948 const ValueT tempX = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0.0, 0.0));
949 const ValueT tempY = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.0, 0.5, 0.0));
950 const ValueT tempZ = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.0, 0.0, 0.5));
951
952 return ValueT(tempX.x(), tempY.y(), tempZ.z());
953}
954
955
956//////////////////////////////////////// StaggeredQuadraticSampler
957
958
959template<class TreeT>
960inline bool
961StaggeredQuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
962 typename TreeT::ValueType& result)
963{
964 using ValueType = typename TreeT::ValueType;
965
966 ValueType tempX, tempY, tempZ;
967 bool active = false;
968
969 active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
970 active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
971 active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
972
973 result.x() = tempX.x();
974 result.y() = tempY.y();
975 result.z() = tempZ.z();
976
977 return active;
978}
979
980template<class TreeT>
981inline typename TreeT::ValueType
982StaggeredQuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
983{
984 using ValueT = typename TreeT::ValueType;
985
986 const ValueT tempX = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0.0, 0.0));
987 const ValueT tempY = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.0, 0.5, 0.0));
988 const ValueT tempZ = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.0, 0.0, 0.5));
989
990 return ValueT(tempX.x(), tempY.y(), tempZ.z());
991}
992
993//////////////////////////////////////// Sampler
994
995template <>
996struct Sampler<0, false> : public PointSampler {};
997
998template <>
999struct Sampler<1, false> : public BoxSampler {};
1000
1001template <>
1002struct Sampler<2, false> : public QuadraticSampler {};
1003
1004template <>
1005struct Sampler<0, true> : public StaggeredPointSampler {};
1006
1007template <>
1008struct Sampler<1, true> : public StaggeredBoxSampler {};
1009
1010template <>
1011struct Sampler<2, true> : public StaggeredQuadraticSampler {};
1012
1013} // namespace tools
1014} // namespace OPENVDB_VERSION_NAME
1015} // namespace openvdb
1016
1017#endif // OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
ValueT value
Definition GridBuilder.h:1290
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition Platform.h:204
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition Platform.h:205
Container class that associates a tree with a transform and metadata.
Definition Grid.h:571
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
Int32 ValueType
Definition Coord.h:32
Definition Transform.h:40
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition Vec3.h:85
T & y()
Definition Vec3.h:86
T & z()
Definition Vec3.h:87
Definition Interpolation.h:544
tools::DualGridSampler< AccT, SamplerT > mSampler
Definition Interpolation.h:574
SamplerT SamlerType
Definition Interpolation.h:550
AlphaMask(const GridT &grid, const MaskT &mask, FloatT min, FloatT max, bool invert)
Definition Interpolation.h:553
FloatT FloatType
Definition Interpolation.h:551
const FloatT mInvNorm
Definition Interpolation.h:575
typename MaskType::ConstAccessor AccT
Definition Interpolation.h:572
const bool mInvert
Definition Interpolation.h:576
MaskT MaskType
Definition Interpolation.h:549
AccT mAcc
Definition Interpolation.h:573
bool operator()(const Coord &xyz, FloatT &a, FloatT &b) const
Definition Interpolation.h:563
GridT GridType
Definition Interpolation.h:548
DualGridSampler(const AccessorType &sourceAccessor, const math::Transform &sourceXform, const math::Transform &targetXform)
ValueAccessor and transform constructor.
Definition Interpolation.h:509
typename TreeT::ValueType ValueType
Definition Interpolation.h:500
bool isAligned() const
Return true if the two grids are aligned.
Definition Interpolation.h:527
typename tree::ValueAccessor< TreeT > AccessorType
Definition Interpolation.h:503
ValueType operator()(const Coord &ijk) const
Return the value of the source grid at the index coordinates, ijk, relative to the target grid.
Definition Interpolation.h:520
This is a simple convenience class that allows for sampling from a source grid into the index space o...
Definition Interpolation.h:446
typename GridOrTreeT::ValueType ValueType
Definition Interpolation.h:448
typename TreeAdapter< GridOrTreeT >::GridType GridType
Definition Interpolation.h:449
bool isAligned() const
Return true if the two grids are aligned.
Definition Interpolation.h:486
DualGridSampler(const GridType &sourceGrid, const math::Transform &targetXform)
Grid and transform constructor.
Definition Interpolation.h:456
DualGridSampler(const TreeType &sourceTree, const math::Transform &sourceXform, const math::Transform &targetXform)
Tree and transform constructor.
Definition Interpolation.h:468
typename TreeAdapter< GridOrTreeT >::TreeType TreeType
Definition Interpolation.h:450
typename TreeAdapter< GridType >::AccessorType AccessorType
Definition Interpolation.h:451
ValueType operator()(const Coord &ijk) const
Return the value of the source grid at the index coordinates, ijk, relative to the target grid (or it...
Definition Interpolation.h:479
typename TreeT::ValueType ValueType
Definition Interpolation.h:369
ValueType isSample(const Coord &ijk) const
Sample value in integer index space.
Definition Interpolation.h:405
typename tree::ValueAccessor< TreeT > AccessorType
Definition Interpolation.h:372
ValueType sampleVoxel(const RealType &x, const RealType &y, const RealType &z) const
Sample a point in index space in the grid.
Definition Interpolation.h:387
ValueType wsSample(const Vec3d &wspoint) const
Sample in world space.
Definition Interpolation.h:418
GridSampler(const AccessorType &acc, const math::Transform &transform)
Definition Interpolation.h:376
SharedPtr< GridSampler > Ptr
Definition Interpolation.h:368
const math::Transform & transform() const
Definition Interpolation.h:380
ValueType sampleVoxel(typename Coord::ValueType i, typename Coord::ValueType j, typename Coord::ValueType k) const
Sample value in integer index space.
Definition Interpolation.h:396
ValueType isSample(const Vec3d &ispoint) const
Sample in fractional index space.
Definition Interpolation.h:409
Class that provides the interface for continuous sampling of values in a tree.
Definition Interpolation.h:284
typename TreeAdapter< GridOrTreeType >::TreeType TreeType
Definition Interpolation.h:289
GridSampler(const GridType &grid)
Definition Interpolation.h:293
typename GridOrTreeType::ValueType ValueType
Definition Interpolation.h:287
ValueType isSample(const Coord &ijk) const
Sample value in integer index space.
Definition Interpolation.h:326
typename TreeAdapter< GridOrTreeType >::GridType GridType
Definition Interpolation.h:288
ValueType sampleVoxel(const RealType &x, const RealType &y, const RealType &z) const
Sample a point in index space in the grid.
Definition Interpolation.h:308
GridSampler(const TreeType &tree, const math::Transform &transform)
Definition Interpolation.h:298
ValueType wsSample(const Vec3d &wspoint) const
Sample in world space.
Definition Interpolation.h:339
typename TreeAdapter< GridOrTreeType >::AccessorType AccessorType
Definition Interpolation.h:290
SharedPtr< GridSampler > Ptr
Definition Interpolation.h:286
const math::Transform & transform() const
Definition Interpolation.h:301
ValueType sampleVoxel(typename Coord::ValueType i, typename Coord::ValueType j, typename Coord::ValueType k) const
Sample value in integer index space.
Definition Interpolation.h:317
ValueType isSample(const Vec3d &ispoint) const
Sample in fractional index space.
Definition Interpolation.h:330
Definition ValueAccessor.h:191
Vec3i ceilVec3(const Vec3R &v)
Definition Interpolation.h:591
Vec3i roundVec3(const Vec3R &v)
Definition Interpolation.h:598
Vec3i floorVec3(const Vec3R &v)
Definition Interpolation.h:584
math::Extrema extrema(const IterT &iter, bool threaded=true)
Iterate over a scalar grid and compute extrema (min/max) of the values of the voxels that are visited...
Definition Statistics.h:354
std::shared_ptr< T > SharedPtr
Definition Types.h:114
Definition Exceptions.h:13
_TreeType TreeType
Definition Grid.h:1061
typename tree::ValueAccessor< TreeType > AccessorType
Definition Grid.h:1072
Definition Interpolation.h:120
static bool staggered()
Definition Interpolation.h:125
static size_t order()
Definition Interpolation.h:126
static int radius()
Definition Interpolation.h:122
static bool mipmap()
Definition Interpolation.h:123
static const char * name()
Definition Interpolation.h:121
static bool consistent()
Definition Interpolation.h:124
Definition Interpolation.h:97
static bool staggered()
Definition Interpolation.h:102
static size_t order()
Definition Interpolation.h:103
static int radius()
Definition Interpolation.h:99
static bool mipmap()
Definition Interpolation.h:100
static const char * name()
Definition Interpolation.h:98
static bool consistent()
Definition Interpolation.h:101
Definition Interpolation.h:163
static bool staggered()
Definition Interpolation.h:168
static size_t order()
Definition Interpolation.h:169
static int radius()
Definition Interpolation.h:165
static bool mipmap()
Definition Interpolation.h:166
static const char * name()
Definition Interpolation.h:164
static bool consistent()
Definition Interpolation.h:167
Provises a unified interface for sampling, i.e. interpolation.
Definition Interpolation.h:64
static const char * name()
static bool sample(const TreeT &inTree, const Vec3R &inCoord, typename TreeT::ValueType &result)
Sample inTree at the floating-point index coordinate inCoord and store the result in result.
static TreeT::ValueType sample(const TreeT &inTree, const Vec3R &inCoord)
Sample inTree at the floating-point index coordinate inCoord.
Definition Interpolation.h:220
static bool staggered()
Definition Interpolation.h:225
static size_t order()
Definition Interpolation.h:226
static int radius()
Definition Interpolation.h:222
static bool mipmap()
Definition Interpolation.h:223
static const char * name()
Definition Interpolation.h:221
static bool consistent()
Definition Interpolation.h:224
Definition Interpolation.h:197
static bool staggered()
Definition Interpolation.h:202
static size_t order()
Definition Interpolation.h:203
static int radius()
Definition Interpolation.h:199
static bool mipmap()
Definition Interpolation.h:200
static const char * name()
Definition Interpolation.h:198
static bool consistent()
Definition Interpolation.h:201
static bool staggered()
Definition Interpolation.h:248
static size_t order()
Definition Interpolation.h:249
static int radius()
Definition Interpolation.h:245
static bool mipmap()
Definition Interpolation.h:246
static const char * name()
Definition Interpolation.h:244
static bool consistent()
Definition Interpolation.h:247
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212