OpenVDB  1.1.0
Interpolation.h
Go to the documentation of this file.
1 
2 //
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 
76 
77 namespace openvdb {
79 namespace OPENVDB_VERSION_NAME {
80 namespace tools {
81 
82 // The following samplers operate in voxel space.
83 // When the samplers are applied to grids holding vector or other non-scalar data,
84 // the data is assumed to be collocated. For example, using the BoxSampler on a grid
85 // with ValueType Vec3f assumes that all three elements in a vector can be assigned
86 // the same physical location.
87 
89 {
90  static const char* name() { return "point"; }
91  static int radius() { return 0; }
92  static bool mipmap() { return false; }
93  static bool consistent() { return true; }
94 
98  template<class TreeT>
99  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
100  typename TreeT::ValueType& result);
101 };
102 
103 
105 {
106  static const char* name() { return "box"; }
107  static int radius() { return 1; }
108  static bool mipmap() { return true; }
109  static bool consistent() { return true; }
110 
114  template<class TreeT>
115  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
116  typename TreeT::ValueType& result);
117 
120  template<class TreeT>
121  static typename TreeT::ValueType sample(const TreeT& inTree, const Vec3R& inCoord);
122 
123 private:
124  template<class ValueT, size_t N>
125  static inline ValueT trilinearInterpolation(ValueT (& data)[N][N][N], const Vec3R& uvw);
126 };
127 
128 
130 {
131  static const char* name() { return "quadratic"; }
132  static int radius() { return 1; }
133  static bool mipmap() { return true; }
134  static bool consistent() { return false; }
135 
139  template<class TreeT>
140  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
141  typename TreeT::ValueType& result);
142 };
143 
144 
146 
147 
148 // The following samplers operate in voxel space and are designed for Vec3
149 // staggered grid data (e.g., fluid simulations using the Marker-and-Cell approach
150 // associate elements of the velocity vector with different physical locations:
151 // the faces of a cube).
152 
154 {
155  static const char* name() { return "point"; }
156  static int radius() { return 0; }
157  static bool mipmap() { return false; }
158  static bool consistent() { return false; }
159 
163  template<class TreeT>
164  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
165  typename TreeT::ValueType& result);
166 };
167 
168 
170 {
171  static const char* name() { return "box"; }
172  static int radius() { return 1; }
173  static bool mipmap() { return true; }
174  static bool consistent() { return false; }
175 
179  template<class TreeT>
180  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
181  typename TreeT::ValueType& result);
182 };
183 
184 
186 {
187  static const char* name() { return "quadratic"; }
188  static int radius() { return 1; }
189  static bool mipmap() { return true; }
190  static bool consistent() { return false; }
191 
195  template<class TreeT>
196  static bool sample(const TreeT& inTree, const Vec3R& inCoord,
197  typename TreeT::ValueType& result);
198 };
199 
200 
202 
203 
208 template<typename TreeOrAccessorType, typename SamplerType>
210 {
211 public:
212  typedef boost::shared_ptr<GridSampler> Ptr;
213  typedef typename TreeOrAccessorType::ValueType ValueType;
214 
218  explicit GridSampler(const TreeOrAccessorType& tree,
219  const math::Transform& transform = math::Transform()):
220  mTree(&tree), mTransform(transform) {}
221 
223 
228  template<typename RealType>
229  ValueType sampleVoxel(const RealType& x, const RealType& y, const RealType& z) const
230  {
231  return isSample(Vec3d(x,y,z));
232  }
233 
236  ValueType isSample(const Vec3d& ispoint) const
237  {
238  ValueType result = zeroVal<ValueType>();
239  SamplerType::sample(*mTree, ispoint, result);
240  return result;
241  }
242 
245  ValueType wsSample(const Vec3d& wspoint) const
246  {
247  ValueType result = zeroVal<ValueType>();
248  SamplerType::sample(*mTree, mTransform.worldToIndex(wspoint), result);
249  return result;
250  }
251 
252 private:
253  const TreeOrAccessorType* mTree;
254  const math::Transform mTransform;
255 };
256 
257 
259 
260 
261 namespace local_util {
262 
263 inline Vec3i
264 floorVec3(const Vec3R& v)
265 {
266  return Vec3i(int(std::floor(v(0))), int(std::floor(v(1))), int(std::floor(v(2))));
267 }
268 
269 
270 inline Vec3i
271 ceilVec3(const Vec3R& v)
272 {
273  return Vec3i(int(std::ceil(v(0))), int(std::ceil(v(1))), int(std::ceil(v(2))));
274 }
275 
276 
277 inline Vec3i
278 roundVec3(const Vec3R& v)
279 {
280  return Vec3i(int(::round(v(0))), int(::round(v(1))), int(::round(v(2))));
281 }
282 
283 } // namespace local_util
284 
285 
287 
288 
289 template<class TreeT>
290 inline bool
291 PointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
292  typename TreeT::ValueType& result)
293 {
294  Vec3i inIdx = local_util::roundVec3(inCoord);
295  return inTree.probeValue(Coord(inIdx), result);
296 }
297 
298 
300 
301 
302 template<class ValueT, size_t N>
303 inline ValueT
304 BoxSampler::trilinearInterpolation(ValueT (& data)[N][N][N], const Vec3R& uvw)
305 {
306  // Trilinear interpolation:
307  // The eight surrounding latice values are used to construct the result. \n
308  // result(x,y,z) =
309  // v000 (1-x)(1-y)(1-z) + v001 (1-x)(1-y)z + v010 (1-x)y(1-z) + v011 (1-x)yz
310  // + v100 x(1-y)(1-z) + v101 x(1-y)z + v110 xy(1-z) + v111 xyz
311 
312  ValueT resultA, resultB;
313 
314  resultA = data[0][0][0] + ValueT((data[0][0][1] - data[0][0][0]) * uvw[2]);
315  resultB = data[0][1][0] + ValueT((data[0][1][1] - data[0][1][0]) * uvw[2]);
316  ValueT result1 = resultA + ValueT((resultB-resultA) * uvw[1]);
317 
318  resultA = data[1][0][0] + ValueT((data[1][0][1] - data[1][0][0]) * uvw[2]);
319  resultB = data[1][1][0] + ValueT((data[1][1][1] - data[1][1][0]) * uvw[2]);
320  ValueT result2 = resultA + ValueT((resultB - resultA) * uvw[1]);
321 
322  return result1 + ValueT(uvw[0] * (result2 - result1));
323 }
324 
325 
326 template<class TreeT>
327 inline bool
328 BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
329  typename TreeT::ValueType& result)
330 {
331  typedef typename TreeT::ValueType ValueT;
332 
333  Vec3i inIdx = local_util::floorVec3(inCoord);
334  Vec3R uvw = inCoord - inIdx;
335 
336  // Retrieve the values of the eight voxels surrounding the
337  // fractional source coordinates.
338  ValueT data[2][2][2];
339 
340  bool hasActiveValues = false;
341  Coord ijk(inIdx);
342  hasActiveValues |= inTree.probeValue(ijk, data[0][0][0]); // i, j, k
343  ijk[2] += 1;
344  hasActiveValues |= inTree.probeValue(ijk, data[0][0][1]); // i, j, k + 1
345  ijk[1] += 1;
346  hasActiveValues |= inTree.probeValue(ijk, data[0][1][1]); // i, j+1, k + 1
347  ijk[2] = inIdx[2];
348  hasActiveValues |= inTree.probeValue(ijk, data[0][1][0]); // i, j+1, k
349  ijk[0] += 1;
350  ijk[1] = inIdx[1];
351  hasActiveValues |= inTree.probeValue(ijk, data[1][0][0]); // i+1, j, k
352  ijk[2] += 1;
353  hasActiveValues |= inTree.probeValue(ijk, data[1][0][1]); // i+1, j, k + 1
354  ijk[1] += 1;
355  hasActiveValues |= inTree.probeValue(ijk, data[1][1][1]); // i+1, j+1, k + 1
356  ijk[2] = inIdx[2];
357  hasActiveValues |= inTree.probeValue(ijk, data[1][1][0]); // i+1, j+1, k
358 
359  result = trilinearInterpolation(data, uvw);
360  return hasActiveValues;
361 }
362 
363 
364 template<class TreeT>
365 inline typename TreeT::ValueType
366 BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord)
367 {
368  typedef typename TreeT::ValueType ValueT;
369 
370  Vec3i inIdx = local_util::floorVec3(inCoord);
371  Vec3R uvw = inCoord - inIdx;
372 
373  // Retrieve the values of the eight voxels surrounding the
374  // fractional source coordinates.
375  ValueT data[2][2][2];
376 
377  Coord ijk(inIdx);
378  data[0][0][0] = inTree.getValue(ijk); // i, j, k
379  ijk[2] += 1;
380  data[0][0][1] = inTree.getValue(ijk); // i, j, k + 1
381  ijk[1] += 1;
382  data[0][1][1] = inTree.getValue(ijk); // i, j+1, k + 1
383  ijk[2] = inIdx[2];
384  data[0][1][0] = inTree.getValue(ijk); // i, j+1, k
385  ijk[0] += 1;
386  ijk[1] = inIdx[1];
387  data[1][0][0] = inTree.getValue(ijk); // i+1, j, k
388  ijk[2] += 1;
389  data[1][0][1] = inTree.getValue(ijk); // i+1, j, k + 1
390  ijk[1] += 1;
391  data[1][1][1] = inTree.getValue(ijk); // i+1, j+1, k + 1
392  ijk[2] = inIdx[2];
393  data[1][1][0] = inTree.getValue(ijk); // i+1, j+1, k
394 
395  return trilinearInterpolation(data, uvw);
396 }
397 
398 
400 
401 
402 template<class TreeT>
403 inline bool
404 QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
405  typename TreeT::ValueType& result)
406 {
407  typedef typename TreeT::ValueType ValueT;
408 
409  Vec3i
410  inIdx = local_util::floorVec3(inCoord),
411  inLoIdx = inIdx - Vec3i(1, 1, 1);
412  Vec3R frac = inCoord - inIdx;
413 
414  // Retrieve the values of the 27 voxels surrounding the
415  // fractional source coordinates.
416  bool active = false;
417  ValueT v[3][3][3];
418  for (int dx = 0, ix = inLoIdx.x(); dx < 3; ++dx, ++ix) {
419  for (int dy = 0, iy = inLoIdx.y(); dy < 3; ++dy, ++iy) {
420  for (int dz = 0, iz = inLoIdx.z(); dz < 3; ++dz, ++iz) {
421  if (inTree.probeValue(Coord(ix, iy, iz), v[dx][dy][dz])) {
422  active = true;
423  }
424  }
425  }
426  }
427 
429  ValueT vx[3];
430  for (int dx = 0; dx < 3; ++dx) {
431  ValueT vy[3];
432  for (int dy = 0; dy < 3; ++dy) {
433  // Fit a parabola to three contiguous samples in z
434  // (at z=-1, z=0 and z=1), then evaluate the parabola at z',
435  // where z' is the fractional part of inCoord.z, i.e.,
436  // inCoord.z - inIdx.z. The coefficients come from solving
437  //
438  // | (-1)^2 -1 1 || a | | v0 |
439  // | 0 0 1 || b | = | v1 |
440  // | 1^2 1 1 || c | | v2 |
441  //
442  // for a, b and c.
443  const ValueT* vz = &v[dx][dy][0];
444  const ValueT
445  az = static_cast<ValueT>(0.5 * (vz[0] + vz[2]) - vz[1]),
446  bz = static_cast<ValueT>(0.5 * (vz[2] - vz[0])),
447  cz = static_cast<ValueT>(vz[1]);
448  vy[dy] = static_cast<ValueT>(frac.z() * (frac.z() * az + bz) + cz);
449  }
450  // Fit a parabola to three interpolated samples in y, then
451  // evaluate the parabola at y', where y' is the fractional
452  // part of inCoord.y.
453  const ValueT
454  ay = static_cast<ValueT>(0.5 * (vy[0] + vy[2]) - vy[1]),
455  by = static_cast<ValueT>(0.5 * (vy[2] - vy[0])),
456  cy = static_cast<ValueT>(vy[1]);
457  vx[dx] = static_cast<ValueT>(frac.y() * (frac.y() * ay + by) + cy);
458  }
459  // Fit a parabola to three interpolated samples in x, then
460  // evaluate the parabola at the fractional part of inCoord.x.
461  const ValueT
462  ax = static_cast<ValueT>(0.5 * (vx[0] + vx[2]) - vx[1]),
463  bx = static_cast<ValueT>(0.5 * (vx[2] - vx[0])),
464  cx = static_cast<ValueT>(vx[1]);
465  result = static_cast<ValueT>(frac.x() * (frac.x() * ax + bx) + cx);
466 
467  return active;
468 }
469 
470 
472 
473 
474 template<class TreeT>
475 inline bool
476 StaggeredPointSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
477  typename TreeT::ValueType& result)
478 {
479  typedef typename TreeT::ValueType ValueType;
480 
481  ValueType tempX, tempY, tempZ;
482  bool active = false;
483 
484  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
485  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
486  active = PointSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
487 
488  result.x() = tempX.x();
489  result.y() = tempY.y();
490  result.z() = tempZ.z();
491 
492  return active;
493 }
494 
495 
497 
498 
499 template<class TreeT>
500 inline bool
501 StaggeredBoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
502  typename TreeT::ValueType& result)
503 {
504  typedef typename TreeT::ValueType ValueType;
505 
506  ValueType tempX, tempY, tempZ;
507  tempX = tempY = tempZ = zeroVal<ValueType>();
508  bool active = false;
509 
510  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
511  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
512  active = BoxSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
513 
514  result.x() = tempX.x();
515  result.y() = tempY.y();
516  result.z() = tempZ.z();
517 
518  return active;
519 }
520 
521 
523 
524 
525 template<class TreeT>
526 inline bool
527 StaggeredQuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
528  typename TreeT::ValueType& result)
529 {
530  typedef typename TreeT::ValueType ValueType;
531 
532  ValueType tempX, tempY, tempZ;
533  bool active = false;
534 
535  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0.5, 0, 0), tempX) || active;
536  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0.5, 0), tempY) || active;
537  active = QuadraticSampler::sample<TreeT>(inTree, inCoord + Vec3R(0, 0, 0.5), tempZ) || active;
538 
539  result.x() = tempX.x();
540  result.y() = tempY.y();
541  result.z() = tempZ.z();
542 
543  return active;
544 }
545 
546 } // namespace tools
547 } // namespace OPENVDB_VERSION_NAME
548 } // namespace openvdb
549 
550 #endif // OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
551 
552 // Copyright (c) 2012-2013 DreamWorks Animation LLC
553 // All rights reserved. This software is distributed under the
554 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )