OpenVDB  0.104.0
Interpolation.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012 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/Platform.h> // for round()
73 
74 
75 namespace openvdb {
77 namespace OPENVDB_VERSION_NAME {
78 namespace tools {
79 
80 // The following samplers operate in voxel space.
81 // When the samplers are applied to grids holding vector or other non-scalar data,
82 // the data is assumed to be collocated. For example, using the BoxSampler on a grid
83 // with ValueType Vec3f assumes that all three elements in a vector can be assigned
84 // the same physical location.
85 
87 {
88  static const char* name() { return "point"; }
89  static int radius() { return 0; }
90  static bool mipmap() { return false; }
91  static bool consistent() { return true; }
92 
96  template<class TreeT>
97  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
98  typename TreeT::ValueType& result);
99 };
100 
101 
103 {
104  static const char* name() { return "box"; }
105  static int radius() { return 1; }
106  static bool mipmap() { return true; }
107  static bool consistent() { return true; }
108 
112  template<class TreeT>
113  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
114  typename TreeT::ValueType& result);
115 
118  template<class TreeT>
119  static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
120 
121 private:
122  template<class ValueT, size_t N>
123  static inline ValueT trilinearInterpolation(ValueT (& data)[N][N][N], const Vec3R& uvw);
124 };
125 
126 
128 {
129  static const char* name() { return "quadratic"; }
130  static int radius() { return 1; }
131  static bool mipmap() { return true; }
132  static bool consistent() { return false; }
133 
137  template<class TreeT>
138  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
139  typename TreeT::ValueType& result);
140 };
141 
142 
144 
145 
146 // The following samplers operate in voxel space and are designed for Vec3
147 // staggered grid data (e.g., fluid simulations using the Marker-and-Cell approach
148 // associate elements of the velocity vector with different physical locations:
149 // the faces of a cube).
150 
152 {
153  static const char* name() { return "point"; }
154  static int radius() { return 0; }
155  static bool mipmap() { return false; }
156  static bool consistent() { return false; }
157 
161  template<class TreeT>
162  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
163  typename TreeT::ValueType& result);
164 };
165 
166 
168 {
169  static const char* name() { return "box"; }
170  static int radius() { return 1; }
171  static bool mipmap() { return true; }
172  static bool consistent() { return false; }
173 
177  template<class TreeT>
178  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
179  typename TreeT::ValueType& result);
180 };
181 
182 
184 {
185  static const char* name() { return "quadratic"; }
186  static int radius() { return 1; }
187  static bool mipmap() { return true; }
188  static bool consistent() { return false; }
189 
193  template<class TreeT>
194  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
195  typename TreeT::ValueType& result);
196 };
197 
198 
200 
201 
206 template<typename TreeOrAccessorType, typename SamplerType>
208 {
209 public:
210  typedef boost::shared_ptr<GridSampler> Ptr;
211  typedef typename TreeOrAccessorType::ValueType ValueType;
212 
216  explicit GridSampler(const TreeOrAccessorType& tree,
217  const math::Transform& transform = math::Transform()):
218  mTree(&tree), mTransform(transform) {}
219 
221 
226  template<typename RealType>
227  ValueType sampleVoxel(const RealType& x, const RealType& y, const RealType& z) const
228  {
229  return isSample(Vec3d(x,y,z));
230  }
231 
234  ValueType isSample(const Vec3d& ispoint) const
235  {
236  ValueType result = zeroVal<ValueType>();
237  SamplerType::sample(*mTree, ispoint, result);
238  return result;
239  }
240 
243  ValueType wsSample(const Vec3d& wspoint) const
244  {
245  ValueType result = zeroVal<ValueType>();
246  SamplerType::sample(*mTree, mTransform.worldToIndex(wspoint), result);
247  return result;
248  }
249 
250 private:
251  const TreeOrAccessorType* mTree;
252  const math::Transform mTransform;
253 };
254 
255 
257 
258 
264 template<typename TreeType>
266 {
267 public:
268  typedef boost::shared_ptr<GridSampling> Ptr;
269  typedef typename TreeType::ValueType ValueType;
271 
273  explicit GridSampling(const TreeType& tree): mTree(&tree) {}
274 
276  virtual ~GridSampling() {}
277 
278  const TreeType& tree() const { return *mTree; }
279 
282  ValueType sampleVoxel(const Vec3R& pt) const { return sampleVoxel(pt[0], pt[1], pt[2]); }
283 
288  virtual ValueType sampleVoxel(Real x, Real y, Real z) const = 0;
289 
290 private:
291  // The grid we are sampling.
292  const TreeType* mTree;
293 };
294 
295 
297 
298 
300 template<typename TreeType>
301 class LinearInterp: public GridSampling<TreeType>
302 {
303 public:
305  typedef boost::shared_ptr<LinearInterp> Ptr;
306  typedef typename TreeType::ValueType ValueType;
307 
309  OPENVDB_DEPRECATED explicit LinearInterp(const TreeType& tree): BaseType(tree) {}
310  virtual ~LinearInterp() {}
311 
317  virtual ValueType sampleVoxel(Real x, Real y, Real z) const
318  {
319  ValueType result = zeroVal<ValueType>();
320  BoxSampler::sample(this->tree(), Vec3R(x, y, z), result);
321  return result;
322  }
323 };
324 
325 
327 
328 
329 template<typename TreeType>
330 class QuadraticInterp: public GridSampling<TreeType>
331 {
332 public:
334  typedef boost::shared_ptr<QuadraticInterp> Ptr;
335  typedef typename TreeType::Ptr TreePtr;
336  typedef typename TreeType::ValueType ValueType;
337 
339  OPENVDB_DEPRECATED explicit QuadraticInterp(const TreeType& tree): BaseType(tree) {}
340  virtual ~QuadraticInterp() {}
341 
343  virtual ValueType sampleVoxel(Real x, Real y, Real z) const
344  {
345  ValueType result = zeroVal<ValueType>();
346  QuadraticSampler::sample(this->tree(), Vec3R(x, y, z), result);
347  return result;
348  }
349 };
350 
351 
353 
354 
355 namespace local_util {
356 
357 inline Vec3i
358 floorVec3(const Vec3R& v)
359 {
360  return Vec3i(int(std::floor(v(0))), int(std::floor(v(1))), int(std::floor(v(2))));
361 }
362 
363 
364 inline Vec3i
365 ceilVec3(const Vec3R& v)
366 {
367  return Vec3i(int(std::ceil(v(0))), int(std::ceil(v(1))), int(std::ceil(v(2))));
368 }
369 
370 
371 inline Vec3i
372 roundVec3(const Vec3R& v)
373 {
374  return Vec3i(int(::round(v(0))), int(::round(v(1))), int(::round(v(2))));
375 }
376 
377 } // namespace local_util
378 
379 
381 
382 
383 template<class TreeT>
384 inline bool
385 PointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
386  typename TreeT::ValueType& result)
387 {
388  Vec3i inIdx = local_util::roundVec3(inCoord);
389  return inTree.probeValue(Coord(inIdx), result);
390 }
391 
392 
394 
395 
396 template<class ValueT, size_t N>
397 inline ValueT
398 BoxSampler::trilinearInterpolation(ValueT (& data)[N][N][N], const Vec3R& uvw)
399 {
400  // Trilinear interpolation:
401  // The eight surrounding latice values are used to construct the result. \n
402  // result(x,y,z) =
403  // v000 (1-x)(1-y)(1-z) + v001 (1-x)(1-y)z + v010 (1-x)y(1-z) + v011 (1-x)yz
404  // + v100 x(1-y)(1-z) + v101 x(1-y)z + v110 xy(1-z) + v111 xyz
405 
406  ValueT resultA, resultB;
407 
408  resultA = data[0][0][0] + ValueT((data[0][0][1] - data[0][0][0]) * uvw[2]);
409  resultB = data[0][1][0] + ValueT((data[0][1][1] - data[0][1][0]) * uvw[2]);
410  ValueT result1 = resultA + ValueT((resultB-resultA) * uvw[1]);
411 
412  resultA = data[1][0][0] + ValueT((data[1][0][1] - data[1][0][0]) * uvw[2]);
413  resultB = data[1][1][0] + ValueT((data[1][1][1] - data[1][1][0]) * uvw[2]);
414  ValueT result2 = resultA + ValueT((resultB - resultA) * uvw[1]);
415 
416  return result1 + ValueT(uvw[0] * (result2 - result1));
417 }
418 
419 
420 template<class TreeT>
421 inline bool
422 BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
423  typename TreeT::ValueType& result)
424 {
425  typedef typename TreeT::ValueType ValueT;
426 
427  Vec3i inIdx = local_util::floorVec3(inCoord);
428  Vec3R uvw = inCoord - inIdx;
429 
430  // Retrieve the values of the eight voxels surrounding the
431  // fractional source coordinates.
432  ValueT data[2][2][2];
433 
434  bool hasActiveValues = false;
435  Coord ijk(inIdx);
436  hasActiveValues |= inTree.probeValue(ijk, data[0][0][0]); // i, j, k
437  ijk[2] += 1;
438  hasActiveValues |= inTree.probeValue(ijk, data[0][0][1]); // i, j, k + 1
439  ijk[1] += 1;
440  hasActiveValues |= inTree.probeValue(ijk, data[0][1][1]); // i, j+1, k + 1
441  ijk[2] = inIdx[2];
442  hasActiveValues |= inTree.probeValue(ijk, data[0][1][0]); // i, j+1, k
443  ijk[0] += 1;
444  ijk[1] = inIdx[1];
445  hasActiveValues |= inTree.probeValue(ijk, data[1][0][0]); // i+1, j, k
446  ijk[2] += 1;
447  hasActiveValues |= inTree.probeValue(ijk, data[1][0][1]); // i+1, j, k + 1
448  ijk[1] += 1;
449  hasActiveValues |= inTree.probeValue(ijk, data[1][1][1]); // i+1, j+1, k + 1
450  ijk[2] = inIdx[2];
451  hasActiveValues |= inTree.probeValue(ijk, data[1][1][0]); // i+1, j+1, k
452 
453  result = trilinearInterpolation(data, uvw);
454  return hasActiveValues;
455 }
456 
457 
458 template<class TreeT>
459 inline typename TreeT::ValueType
460 BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
461 {
462  typedef typename TreeT::ValueType ValueT;
463 
464  Vec3i inIdx = local_util::floorVec3(inCoord);
465  Vec3R uvw = inCoord - inIdx;
466 
467  // Retrieve the values of the eight voxels surrounding the
468  // fractional source coordinates.
469  ValueT data[2][2][2];
470 
471  Coord ijk(inIdx);
472  data[0][0][0] = inTree.getValue(ijk); // i, j, k
473  ijk[2] += 1;
474  data[0][0][1] = inTree.getValue(ijk); // i, j, k + 1
475  ijk[1] += 1;
476  data[0][1][1] = inTree.getValue(ijk); // i, j+1, k + 1
477  ijk[2] = inIdx[2];
478  data[0][1][0] = inTree.getValue(ijk); // i, j+1, k
479  ijk[0] += 1;
480  ijk[1] = inIdx[1];
481  data[1][0][0] = inTree.getValue(ijk); // i+1, j, k
482  ijk[2] += 1;
483  data[1][0][1] = inTree.getValue(ijk); // i+1, j, k + 1
484  ijk[1] += 1;
485  data[1][1][1] = inTree.getValue(ijk); // i+1, j+1, k + 1
486  ijk[2] = inIdx[2];
487  data[1][1][0] = inTree.getValue(ijk); // i+1, j+1, k
488 
489  return trilinearInterpolation(data, uvw);
490 }
491 
492 
494 
495 
496 template<class TreeT>
497 inline bool
498 QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
499  typename TreeT::ValueType& result)
500 {
501  typedef typename TreeT::ValueType ValueT;
502 
503  Vec3i
504  inIdx = local_util::floorVec3(inCoord),
505  inLoIdx = inIdx - Vec3i(1, 1, 1);
506  Vec3R frac = inCoord - inIdx;
507 
508  // Retrieve the values of the 27 voxels surrounding the
509  // fractional source coordinates.
510  bool active = false;
511  ValueT v[3][3][3];
512  for (int dx = 0, ix = inLoIdx.x(); dx < 3; ++dx, ++ix) {
513  for (int dy = 0, iy = inLoIdx.y(); dy < 3; ++dy, ++iy) {
514  for (int dz = 0, iz = inLoIdx.z(); dz < 3; ++dz, ++iz) {
515  if (inTree.probeValue(Coord(ix, iy, iz), v[dx][dy][dz])) {
516  active = true;
517  }
518  }
519  }
520  }
521 
523  ValueT vx[3];
524  for (int dx = 0; dx < 3; ++dx) {
525  ValueT vy[3];
526  for (int dy = 0; dy < 3; ++dy) {
527  // Fit a parabola to three contiguous samples in z
528  // (at z=-1, z=0 and z=1), then evaluate the parabola at z',
529  // where z' is the fractional part of inCoord.z, i.e.,
530  // inCoord.z - inIdx.z. The coefficients come from solving
531  //
532  // | (-1)^2 -1 1 || a | | v0 |
533  // | 0 0 1 || b | = | v1 |
534  // | 1^2 1 1 || c | | v2 |
535  //
536  // for a, b and c.
537  const ValueT* vz = &v[dx][dy][0];
538  const ValueT
539  az = static_cast<ValueT>(0.5 * (vz[0] + vz[2]) - vz[1]),
540  bz = static_cast<ValueT>(0.5 * (vz[2] - vz[0])),
541  cz = static_cast<ValueT>(vz[1]);
542  vy[dy] = static_cast<ValueT>(frac.z() * (frac.z() * az + bz) + cz);
543  }
544  // Fit a parabola to three interpolated samples in y, then
545  // evaluate the parabola at y', where y' is the fractional
546  // part of inCoord.y.
547  const ValueT
548  ay = static_cast<ValueT>(0.5 * (vy[0] + vy[2]) - vy[1]),
549  by = static_cast<ValueT>(0.5 * (vy[2] - vy[0])),
550  cy = static_cast<ValueT>(vy[1]);
551  vx[dx] = static_cast<ValueT>(frac.y() * (frac.y() * ay + by) + cy);
552  }
553  // Fit a parabola to three interpolated samples in x, then
554  // evaluate the parabola at the fractional part of inCoord.x.
555  const ValueT
556  ax = static_cast<ValueT>(0.5 * (vx[0] + vx[2]) - vx[1]),
557  bx = static_cast<ValueT>(0.5 * (vx[2] - vx[0])),
558  cx = static_cast<ValueT>(vx[1]);
559  result = static_cast<ValueT>(frac.x() * (frac.x() * ax + bx) + cx);
560 
561  return active;
562 }
563 
564 
566 
567 
568 template<class TreeT>
569 inline bool
570 StaggeredPointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
571  typename TreeT::ValueType& result)
572 {
573  typedef typename TreeT::ValueType ValueType;
574 
575  ValueType tempX, tempY, tempZ;
576  bool active = false;
577 
578  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
579  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
580  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
581 
582  result.x() = tempX.x();
583  result.y() = tempY.y();
584  result.z() = tempZ.z();
585 
586  return active;
587 }
588 
589 
591 
592 
593 template<class TreeT>
594 inline bool
595 StaggeredBoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
596  typename TreeT::ValueType& result)
597 {
598  typedef typename TreeT::ValueType ValueType;
599 
600  ValueType tempX, tempY, tempZ;
601  tempX = tempY = tempZ = zeroVal<ValueType>();
602  bool active = false;
603 
604  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
605  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
606  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
607 
608  result.x() = tempX.x();
609  result.y() = tempY.y();
610  result.z() = tempZ.z();
611 
612  return active;
613 }
614 
615 
617 
618 
619 template<class TreeT>
620 inline bool
621 StaggeredQuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
622  typename TreeT::ValueType& result)
623 {
624  typedef typename TreeT::ValueType ValueType;
625 
626  ValueType tempX, tempY, tempZ;
627  bool active = false;
628 
629  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
630  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
631  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
632 
633  result.x() = tempX.x();
634  result.y() = tempY.y();
635  result.z() = tempZ.z();
636 
637  return active;
638 }
639 
640 } // namespace tools
641 } // namespace OPENVDB_VERSION_NAME
642 } // namespace openvdb
643 
644 #endif // OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
645 
646 // Copyright (c) 2012 DreamWorks Animation LLC
647 // All rights reserved. This software is distributed under the
648 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )