OpenVDB  9.0.1
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 
52 namespace openvdb {
54 namespace OPENVDB_VERSION_NAME {
55 namespace 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
62 template <size_t Order, bool Staggered = false>
63 struct 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).
282 template<typename GridOrTreeType, typename SamplerType>
284 {
285 public:
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
317  ValueType sampleVoxel(typename Coord::ValueType i,
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 
346 private:
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.
364 template<typename TreeT, typename SamplerType>
365 class GridSampler<tree::ValueAccessor<TreeT>, SamplerType>
366 {
367 public:
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
396  ValueType sampleVoxel(typename Coord::ValueType i,
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 
425 private:
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!
443 template<typename GridOrTreeT,
444  typename SamplerT>
446 {
447 public:
448  using ValueType = typename GridOrTreeT::ValueType;
450  using TreeType = typename TreeAdapter<GridOrTreeT>::TreeType;
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; }
487 private:
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.
495 template<typename TreeT,
496  typename SamplerT>
497 class 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; }
528 private:
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
539 template <typename GridT,
540  typename MaskT,
541  typename SamplerT = tools::BoxSampler,
542  typename FloatT = float>
544 {
545 public:
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 
571 protected:
572  using AccT = typename MaskType::ConstAccessor;
575  const FloatT mMin, mInvNorm;
576  const bool mInvert;
577 };// AlphaMask
578 
579 ////////////////////////////////////////
580 
581 namespace local_util {
582 
583 inline Vec3i
584 floorVec3(const Vec3R& v)
585 {
586  return Vec3i(int(std::floor(v(0))), int(std::floor(v(1))), int(std::floor(v(2))));
587 }
588 
589 
590 inline Vec3i
591 ceilVec3(const Vec3R& v)
592 {
593  return Vec3i(int(std::ceil(v(0))), int(std::ceil(v(1))), int(std::ceil(v(2))));
594 }
595 
596 
597 inline Vec3i
598 roundVec3(const Vec3R& v)
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 
609 template<class TreeT>
610 inline bool
611 PointSampler::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 
617 template<class TreeT>
618 inline typename TreeT::ValueType
619 PointSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
620 {
621  return inTree.getValue(Coord(local_util::roundVec3(inCoord)));
622 }
623 
624 
625 //////////////////////////////////////// BoxSampler
626 
627 template<class ValueT, class TreeT, size_t N>
628 inline void
629 BoxSampler::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 
656 template<class ValueT, class TreeT, size_t N>
657 inline bool
658 BoxSampler::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 
688 template<class ValueT, size_t N>
689 inline void
690 BoxSampler::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 
710 template<class ValueT, size_t N>
711 inline ValueT
712 BoxSampler::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 
741 template<class TreeT>
742 inline bool
743 BoxSampler::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 
763 template<class TreeT>
764 inline typename TreeT::ValueType
765 BoxSampler::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 
784 template<class ValueT, size_t N>
785 inline ValueT
786 QuadraticSampler::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 
828 template<class TreeT>
829 inline bool
830 QuadraticSampler::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 
855 template<class TreeT>
856 inline typename TreeT::ValueType
857 QuadraticSampler::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 
882 template<class TreeT>
883 inline bool
884 StaggeredPointSampler::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 
903 template<class TreeT>
904 inline typename TreeT::ValueType
905 StaggeredPointSampler::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 
920 template<class TreeT>
921 inline bool
922 StaggeredBoxSampler::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 
942 template<class TreeT>
943 inline typename TreeT::ValueType
944 StaggeredBoxSampler::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 
959 template<class TreeT>
960 inline bool
961 StaggeredQuadraticSampler::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 
980 template<class TreeT>
981 inline typename TreeT::ValueType
982 StaggeredQuadraticSampler::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 
995 template <>
996 struct Sampler<0, false> : public PointSampler {};
997 
998 template <>
999 struct Sampler<1, false> : public BoxSampler {};
1000 
1001 template <>
1002 struct Sampler<2, false> : public QuadraticSampler {};
1003 
1004 template <>
1005 struct Sampler<0, true> : public StaggeredPointSampler {};
1006 
1007 template <>
1008 struct Sampler<1, true> : public StaggeredBoxSampler {};
1009 
1010 template <>
1011 struct 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
Definition: Interpolation.h:196
#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:206
GridSampler(const TreeType &tree, const math::Transform &transform)
Definition: Interpolation.h:298
static bool staggered()
Definition: Interpolation.h:225
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x².
Definition: Math.h:286
Vec3i ceilVec3(const Vec3R &v)
Definition: Interpolation.h:591
Class that provides the interface for continuous sampling of values in a tree.
Definition: Interpolation.h:283
static bool mipmap()
Definition: Interpolation.h:166
static bool mipmap()
Definition: Interpolation.h:123
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
static bool staggered()
Definition: Interpolation.h:168
static bool consistent()
Definition: Interpolation.h:167
MaskT MaskType
Definition: Interpolation.h:549
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
static int radius()
Definition: Interpolation.h:199
static bool mipmap()
Definition: Interpolation.h:100
static bool staggered()
Definition: Interpolation.h:102
Definition: Interpolation.h:543
SamplerT SamlerType
Definition: Interpolation.h:550
static int radius()
Definition: Interpolation.h:99
Vec3< int32_t > Vec3i
Definition: Vec3.h:665
T & z()
Definition: Vec3.h:91
static int radius()
Definition: Interpolation.h:222
T & y()
Definition: Vec3.h:90
typename tree::ValueAccessor< TreeT > AccessorType
Definition: Interpolation.h:372
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
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
SharedPtr< GridSampler > Ptr
Definition: Interpolation.h:368
static const char * name()
Definition: Interpolation.h:198
static int radius()
Definition: Interpolation.h:122
GridSampler(const GridType &grid)
Definition: Interpolation.h:293
typename GridOrTreeType::ValueType ValueType
Definition: Interpolation.h:287
Definition: Transform.h:39
std::shared_ptr< T > SharedPtr
Definition: Types.h:114
static size_t order()
Definition: Interpolation.h:226
const FloatT mMin
Definition: Interpolation.h:575
static bool consistent()
Definition: Interpolation.h:101
const math::Transform & transform() const
Definition: Interpolation.h:380
typename tree::ValueAccessor< TreeT > AccessorType
Definition: Interpolation.h:503
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:207
ValueType isSample(const Coord &ijk) const
Sample value in integer index space.
Definition: Interpolation.h:405
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
static bool mipmap()
Definition: Interpolation.h:200
static size_t order()
Definition: Interpolation.h:203
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
static const char * name()
Definition: Interpolation.h:221
tools::DualGridSampler< AccT, SamplerT > mSampler
Definition: Interpolation.h:574
typename TreeT::ValueType ValueType
Definition: Interpolation.h:500
static size_t order()
Definition: Interpolation.h:126
static int radius()
Definition: Interpolation.h:165
ValueType wsSample(const Vec3d &wspoint) const
Sample in world space.
Definition: Interpolation.h:339
DualGridSampler(const GridType &sourceGrid, const math::Transform &targetXform)
Grid and transform constructor.
Definition: Interpolation.h:456
typename TreeT::ValueType ValueType
Definition: Interpolation.h:369
static bool consistent()
Definition: Interpolation.h:247
bool isAligned() const
Return true if the two grids are aligned.
Definition: Interpolation.h:486
static size_t order()
Definition: Interpolation.h:169
This is a simple convenience class that allows for sampling from a source grid into the index space o...
Definition: Interpolation.h:445
math::Vec3< Real > Vec3R
Definition: Types.h:72
typename TreeAdapter< GridOrTreeType >::TreeType TreeType
Definition: Interpolation.h:289
GridSampler(const AccessorType &acc, const math::Transform &transform)
Definition: Interpolation.h:376
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
DualGridSampler(const TreeType &sourceTree, const math::Transform &sourceXform, const math::Transform &targetXform)
Tree and transform constructor.
Definition: Interpolation.h:468
static bool staggered()
Definition: Interpolation.h:248
Definition: Exceptions.h:13
const bool mInvert
Definition: Interpolation.h:576
ValueT value
Definition: GridBuilder.h:1287
Definition: Interpolation.h:119
bool operator()(const Coord &xyz, FloatT &a, FloatT &b) const
Definition: Interpolation.h:563
static int radius()
Definition: Interpolation.h:245
typename TreeAdapter< GridOrTreeType >::GridType GridType
Definition: Interpolation.h:288
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:598
static size_t order()
Definition: Interpolation.h:103
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:659
const math::Transform & transform() const
Definition: Interpolation.h:301
static bool staggered()
Definition: Interpolation.h:125
Vec3i roundVec3(const Vec3R &v)
Definition: Interpolation.h:598
Vec3i floorVec3(const Vec3R &v)
Definition: Interpolation.h:584
AccT mAcc
Definition: Interpolation.h:573
static bool consistent()
Definition: Interpolation.h:224
bool isAligned() const
Return true if the two grids are aligned.
Definition: Interpolation.h:527
Definition: Interpolation.h:219
ValueType wsSample(const Vec3d &wspoint) const
Sample in world space.
Definition: Interpolation.h:418
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
ValueType isSample(const Vec3d &ispoint) const
Sample in fractional index space.
Definition: Interpolation.h:330
static const char * name()
Definition: Interpolation.h:121
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:89
GridT GridType
Definition: Interpolation.h:548
Definition: Interpolation.h:96
static bool staggered()
Definition: Interpolation.h:202
static bool consistent()
Definition: Interpolation.h:201
Provises a unified interface for sampling, i.e. interpolation.
Definition: Interpolation.h:63
_TreeType TreeType
Definition: Grid.h:1072
typename TreeAdapter< GridType >::AccessorType AccessorType
Definition: Interpolation.h:451
typename TreeAdapter< GridOrTreeType >::AccessorType AccessorType
Definition: Interpolation.h:290
ValueType isSample(const Vec3d &ispoint) const
Sample in fractional index space.
Definition: Interpolation.h:409
FloatT FloatType
Definition: Interpolation.h:551
static bool mipmap()
Definition: Interpolation.h:223
typename TreeAdapter< AccT >::GridType GridType
Definition: Interpolation.h:449
static bool consistent()
Definition: Interpolation.h:124
typename tree::ValueAccessor< TreeType > AccessorType
Definition: Grid.h:1083
static const char * name()
Definition: Interpolation.h:244
typename AccT::ValueType ValueType
Definition: Interpolation.h:448
AlphaMask(const GridT &grid, const MaskT &mask, FloatT min, FloatT max, bool invert)
Definition: Interpolation.h:553
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
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
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
SharedPtr< GridSampler > Ptr
Definition: Interpolation.h:286
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
static const char * name()
Definition: Interpolation.h:98
Definition: Interpolation.h:162
static bool mipmap()
Definition: Interpolation.h:246
DualGridSampler(const AccessorType &sourceAccessor, const math::Transform &sourceXform, const math::Transform &targetXform)
ValueAccessor and transform constructor.
Definition: Interpolation.h:509
static size_t order()
Definition: Interpolation.h:249
ValueType isSample(const Coord &ijk) const
Sample value in integer index space.
Definition: Interpolation.h:326
static const char * name()
Definition: Interpolation.h:164
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202
typename MaskType::ConstAccessor AccT
Definition: Interpolation.h:572