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