OpenVDB  1.2.0
Interpolation.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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 //
30 //
66 
67 #ifndef OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
68 #define OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
69 
70 #include <cmath>
71 #include <boost/shared_ptr.hpp>
72 #include <openvdb/version.h> // for OPENVDB_VERSION_NAME
73 #include <openvdb/Platform.h> // for round()
74 #include <openvdb/math/Transform.h> // for Transform
75 #include <openvdb/Grid.h>
76 #include <openvdb/tree/ValueAccessor.h>
77 
78 namespace openvdb {
80 namespace OPENVDB_VERSION_NAME {
81 namespace tools {
82 
83 // The following samplers operate in voxel space.
84 // When the samplers are applied to grids holding vector or other non-scalar data,
85 // the data is assumed to be collocated. For example, using the BoxSampler on a grid
86 // with ValueType Vec3f assumes that all three elements in a vector can be assigned
87 // the same physical location.
88 
90 {
91  static const char* name() { return "point"; }
92  static int radius() { return 0; }
93  static bool mipmap() { return false; }
94  static bool consistent() { return true; }
95 
99  template<class TreeT>
100  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
101  typename TreeT::ValueType& result);
102 };
103 
104 
106 {
107  static const char* name() { return "box"; }
108  static int radius() { return 1; }
109  static bool mipmap() { return true; }
110  static bool consistent() { return true; }
111 
115  template<class TreeT>
116  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
117  typename TreeT::ValueType& result);
118 
121  template<class TreeT>
122  static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
123 
124 private:
125  template<class ValueT, size_t N>
126  static inline ValueT trilinearInterpolation(ValueT (& data)[N][N][N], const Vec3R& uvw);
127 };
128 
129 
131 {
132  static const char* name() { return "quadratic"; }
133  static int radius() { return 1; }
134  static bool mipmap() { return true; }
135  static bool consistent() { return false; }
136 
140  template<class TreeT>
141  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
142  typename TreeT::ValueType& result);
143 };
144 
145 
147 
148 
149 // The following samplers operate in voxel space and are designed for Vec3
150 // staggered grid data (e.g., fluid simulations using the Marker-and-Cell approach
151 // associate elements of the velocity vector with different physical locations:
152 // the faces of a cube).
153 
155 {
156  static const char* name() { return "point"; }
157  static int radius() { return 0; }
158  static bool mipmap() { return false; }
159  static bool consistent() { return false; }
160 
164  template<class TreeT>
165  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
166  typename TreeT::ValueType& result);
167 };
168 
169 
171 {
172  static const char* name() { return "box"; }
173  static int radius() { return 1; }
174  static bool mipmap() { return true; }
175  static bool consistent() { return false; }
176 
180  template<class TreeT>
181  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
182  typename TreeT::ValueType& result);
183 };
184 
185 
187 {
188  static const char* name() { return "quadratic"; }
189  static int radius() { return 1; }
190  static bool mipmap() { return true; }
191  static bool consistent() { return false; }
192 
196  template<class TreeT>
197  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
198  typename TreeT::ValueType& result);
199 };
200 
201 
203 
204 
219 template<typename GridOrTreeType, typename SamplerType>
221 {
222 public:
223  typedef boost::shared_ptr<GridSampler> Ptr;
224  typedef typename GridOrTreeType::ValueType ValueType;
228 
230  explicit GridSampler(const GridType& grid)
231  : mTree(&(grid.tree())), mTransform(&(grid.transform())) {}
232 
235  GridSampler(const TreeType& tree, const math::Transform& transform)
236  : mTree(&tree), mTransform(&transform) {}
237 
238  const math::Transform& transform() const { return *mTransform; }
239 
244  template<typename RealType>
245  ValueType sampleVoxel(const RealType& x, const RealType& y, const RealType& z) const
246  {
247  return this->isSample(Vec3d(x,y,z));
248  }
249 
255  typename Coord::ValueType j,
256  typename Coord::ValueType k) const
257  {
258  return this->isSample(Coord(i,j,k));
259  }
260 
263  ValueType isSample(const Coord& ijk) const { return mTree->getValue(ijk); }
264 
267  ValueType isSample(const Vec3d& ispoint) const
268  {
269  ValueType result = zeroVal<ValueType>();
270  SamplerType::sample(*mTree, ispoint, result);
271  return result;
272  }
273 
276  ValueType wsSample(const Vec3d& wspoint) const
277  {
278  ValueType result = zeroVal<ValueType>();
279  SamplerType::sample(*mTree, mTransform->worldToIndex(wspoint), result);
280  return result;
281  }
282 
283 private:
284  const TreeType* mTree;
285  const math::Transform* mTransform;
286 }; // class GridSampler
287 
288 
301 template<typename TreeT, typename SamplerType>
302 class GridSampler<tree::ValueAccessor<TreeT>, SamplerType>
303 {
304 public:
305  typedef boost::shared_ptr<GridSampler> Ptr;
306  typedef typename TreeT::ValueType ValueType;
307  typedef TreeT TreeType;
310 
313  GridSampler(const AccessorType& acc, const math::Transform& transform)
314  : mAccessor(&acc), mTransform(&transform) {}
315 
316  const math::Transform& transform() const { return *mTransform; }
317 
322  template<typename RealType>
323  ValueType sampleVoxel(const RealType& x, const RealType& y, const RealType& z) const
324  {
325  return this->isSample(Vec3d(x,y,z));
326  }
327 
333  typename Coord::ValueType j,
334  typename Coord::ValueType k) const
335  {
336  return this->isSample(Coord(i,j,k));
337  }
338 
341  ValueType isSample(const Coord& ijk) const { return mAccessor->getValue(ijk); }
342 
345  ValueType isSample(const Vec3d& ispoint) const
346  {
347  ValueType result = zeroVal<ValueType>();
348  SamplerType::sample(*mAccessor, ispoint, result);
349  return result;
350  }
351 
354  ValueType wsSample(const Vec3d& wspoint) const
355  {
356  ValueType result = zeroVal<ValueType>();
357  SamplerType::sample(*mAccessor, mTransform->worldToIndex(wspoint), result);
358  return result;
359  }
360 
361 private:
362  const AccessorType* mAccessor;//not thread-safe!
363  const math::Transform* mTransform;
364 };
365 
366 
368 
369 
370 namespace local_util {
371 
372 inline Vec3i
373 floorVec3(const Vec3R& v)
374 {
375  return Vec3i(int(std::floor(v(0))), int(std::floor(v(1))), int(std::floor(v(2))));
376 }
377 
378 
379 inline Vec3i
380 ceilVec3(const Vec3R& v)
381 {
382  return Vec3i(int(std::ceil(v(0))), int(std::ceil(v(1))), int(std::ceil(v(2))));
383 }
384 
385 
386 inline Vec3i
387 roundVec3(const Vec3R& v)
388 {
389  return Vec3i(int(::round(v(0))), int(::round(v(1))), int(::round(v(2))));
390 }
391 
392 } // namespace local_util
393 
394 
396 
397 
398 template<class TreeT>
399 inline bool
400 PointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
401  typename TreeT::ValueType& result)
402 {
403  Vec3i inIdx = local_util::roundVec3(inCoord);
404  return inTree.probeValue(Coord(inIdx), result);
405 }
406 
407 
409 
410 
411 template<class ValueT, size_t N>
412 inline ValueT
413 BoxSampler::trilinearInterpolation(ValueT (& data)[N][N][N], const Vec3R& uvw)
414 {
415  // Trilinear interpolation:
416  // The eight surrounding latice values are used to construct the result. \n
417  // result(x,y,z) =
418  // v000 (1-x)(1-y)(1-z) + v001 (1-x)(1-y)z + v010 (1-x)y(1-z) + v011 (1-x)yz
419  // + v100 x(1-y)(1-z) + v101 x(1-y)z + v110 xy(1-z) + v111 xyz
420 
421  ValueT resultA, resultB;
422 
423  resultA = data[0][0][0] + ValueT((data[0][0][1] - data[0][0][0]) * uvw[2]);
424  resultB = data[0][1][0] + ValueT((data[0][1][1] - data[0][1][0]) * uvw[2]);
425  ValueT result1 = resultA + ValueT((resultB-resultA) * uvw[1]);
426 
427  resultA = data[1][0][0] + ValueT((data[1][0][1] - data[1][0][0]) * uvw[2]);
428  resultB = data[1][1][0] + ValueT((data[1][1][1] - data[1][1][0]) * uvw[2]);
429  ValueT result2 = resultA + ValueT((resultB - resultA) * uvw[1]);
430 
431  return result1 + ValueT(uvw[0] * (result2 - result1));
432 }
433 
434 
435 template<class TreeT>
436 inline bool
437 BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
438  typename TreeT::ValueType& result)
439 {
440  typedef typename TreeT::ValueType ValueT;
441 
442  Vec3i inIdx = local_util::floorVec3(inCoord);
443  Vec3R uvw = inCoord - inIdx;
444 
445  // Retrieve the values of the eight voxels surrounding the
446  // fractional source coordinates.
447  ValueT data[2][2][2];
448 
449  bool hasActiveValues = false;
450  Coord ijk(inIdx);
451  hasActiveValues |= inTree.probeValue(ijk, data[0][0][0]); // i, j, k
452  ijk[2] += 1;
453  hasActiveValues |= inTree.probeValue(ijk, data[0][0][1]); // i, j, k + 1
454  ijk[1] += 1;
455  hasActiveValues |= inTree.probeValue(ijk, data[0][1][1]); // i, j+1, k + 1
456  ijk[2] = inIdx[2];
457  hasActiveValues |= inTree.probeValue(ijk, data[0][1][0]); // i, j+1, k
458  ijk[0] += 1;
459  ijk[1] = inIdx[1];
460  hasActiveValues |= inTree.probeValue(ijk, data[1][0][0]); // i+1, j, k
461  ijk[2] += 1;
462  hasActiveValues |= inTree.probeValue(ijk, data[1][0][1]); // i+1, j, k + 1
463  ijk[1] += 1;
464  hasActiveValues |= inTree.probeValue(ijk, data[1][1][1]); // i+1, j+1, k + 1
465  ijk[2] = inIdx[2];
466  hasActiveValues |= inTree.probeValue(ijk, data[1][1][0]); // i+1, j+1, k
467 
468  result = trilinearInterpolation(data, uvw);
469  return hasActiveValues;
470 }
471 
472 
473 template<class TreeT>
474 inline typename TreeT::ValueType
475 BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
476 {
477  typedef typename TreeT::ValueType ValueT;
478 
479  Vec3i inIdx = local_util::floorVec3(inCoord);
480  Vec3R uvw = inCoord - inIdx;
481 
482  // Retrieve the values of the eight voxels surrounding the
483  // fractional source coordinates.
484  ValueT data[2][2][2];
485 
486  Coord ijk(inIdx);
487  data[0][0][0] = inTree.getValue(ijk); // i, j, k
488  ijk[2] += 1;
489  data[0][0][1] = inTree.getValue(ijk); // i, j, k + 1
490  ijk[1] += 1;
491  data[0][1][1] = inTree.getValue(ijk); // i, j+1, k + 1
492  ijk[2] = inIdx[2];
493  data[0][1][0] = inTree.getValue(ijk); // i, j+1, k
494  ijk[0] += 1;
495  ijk[1] = inIdx[1];
496  data[1][0][0] = inTree.getValue(ijk); // i+1, j, k
497  ijk[2] += 1;
498  data[1][0][1] = inTree.getValue(ijk); // i+1, j, k + 1
499  ijk[1] += 1;
500  data[1][1][1] = inTree.getValue(ijk); // i+1, j+1, k + 1
501  ijk[2] = inIdx[2];
502  data[1][1][0] = inTree.getValue(ijk); // i+1, j+1, k
503 
504  return trilinearInterpolation(data, uvw);
505 }
506 
507 
509 
510 
511 template<class TreeT>
512 inline bool
513 QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
514  typename TreeT::ValueType& result)
515 {
516  typedef typename TreeT::ValueType ValueT;
517 
518  Vec3i
519  inIdx = local_util::floorVec3(inCoord),
520  inLoIdx = inIdx - Vec3i(1, 1, 1);
521  Vec3R frac = inCoord - inIdx;
522 
523  // Retrieve the values of the 27 voxels surrounding the
524  // fractional source coordinates.
525  bool active = false;
526  ValueT v[3][3][3];
527  for (int dx = 0, ix = inLoIdx.x(); dx < 3; ++dx, ++ix) {
528  for (int dy = 0, iy = inLoIdx.y(); dy < 3; ++dy, ++iy) {
529  for (int dz = 0, iz = inLoIdx.z(); dz < 3; ++dz, ++iz) {
530  if (inTree.probeValue(Coord(ix, iy, iz), v[dx][dy][dz])) {
531  active = true;
532  }
533  }
534  }
535  }
536 
538  ValueT vx[3];
539  for (int dx = 0; dx < 3; ++dx) {
540  ValueT vy[3];
541  for (int dy = 0; dy < 3; ++dy) {
542  // Fit a parabola to three contiguous samples in z
543  // (at z=-1, z=0 and z=1), then evaluate the parabola at z',
544  // where z' is the fractional part of inCoord.z, i.e.,
545  // inCoord.z - inIdx.z. The coefficients come from solving
546  //
547  // | (-1)^2 -1 1 || a | | v0 |
548  // | 0 0 1 || b | = | v1 |
549  // | 1^2 1 1 || c | | v2 |
550  //
551  // for a, b and c.
552  const ValueT* vz = &v[dx][dy][0];
553  const ValueT
554  az = static_cast<ValueT>(0.5 * (vz[0] + vz[2]) - vz[1]),
555  bz = static_cast<ValueT>(0.5 * (vz[2] - vz[0])),
556  cz = static_cast<ValueT>(vz[1]);
557  vy[dy] = static_cast<ValueT>(frac.z() * (frac.z() * az + bz) + cz);
558  }
559  // Fit a parabola to three interpolated samples in y, then
560  // evaluate the parabola at y', where y' is the fractional
561  // part of inCoord.y.
562  const ValueT
563  ay = static_cast<ValueT>(0.5 * (vy[0] + vy[2]) - vy[1]),
564  by = static_cast<ValueT>(0.5 * (vy[2] - vy[0])),
565  cy = static_cast<ValueT>(vy[1]);
566  vx[dx] = static_cast<ValueT>(frac.y() * (frac.y() * ay + by) + cy);
567  }
568  // Fit a parabola to three interpolated samples in x, then
569  // evaluate the parabola at the fractional part of inCoord.x.
570  const ValueT
571  ax = static_cast<ValueT>(0.5 * (vx[0] + vx[2]) - vx[1]),
572  bx = static_cast<ValueT>(0.5 * (vx[2] - vx[0])),
573  cx = static_cast<ValueT>(vx[1]);
574  result = static_cast<ValueT>(frac.x() * (frac.x() * ax + bx) + cx);
575 
576  return active;
577 }
578 
579 
581 
582 
583 template<class TreeT>
584 inline bool
585 StaggeredPointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
586  typename TreeT::ValueType& result)
587 {
588  typedef typename TreeT::ValueType ValueType;
589 
590  ValueType tempX, tempY, tempZ;
591  bool active = false;
592 
593  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
594  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
595  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
596 
597  result.x() = tempX.x();
598  result.y() = tempY.y();
599  result.z() = tempZ.z();
600 
601  return active;
602 }
603 
604 
606 
607 
608 template<class TreeT>
609 inline bool
610 StaggeredBoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
611  typename TreeT::ValueType& result)
612 {
613  typedef typename TreeT::ValueType ValueType;
614 
615  ValueType tempX, tempY, tempZ;
616  tempX = tempY = tempZ = zeroVal<ValueType>();
617  bool active = false;
618 
619  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
620  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
621  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
622 
623  result.x() = tempX.x();
624  result.y() = tempY.y();
625  result.z() = tempZ.z();
626 
627  return active;
628 }
629 
630 
632 
633 
634 template<class TreeT>
635 inline bool
636 StaggeredQuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
637  typename TreeT::ValueType& result)
638 {
639  typedef typename TreeT::ValueType ValueType;
640 
641  ValueType tempX, tempY, tempZ;
642  bool active = false;
643 
644  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
645  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
646  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
647 
648  result.x() = tempX.x();
649  result.y() = tempY.y();
650  result.z() = tempZ.z();
651 
652  return active;
653 }
654 
655 } // namespace tools
656 } // namespace OPENVDB_VERSION_NAME
657 } // namespace openvdb
658 
659 #endif // OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
660 
661 // Copyright (c) 2012-2013 DreamWorks Animation LLC
662 // All rights reserved. This software is distributed under the
663 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )