OpenVDB  3.0.0
LevelSetAdvect.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2014 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 //
32 //
38 
39 #ifndef OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
40 #define OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
41 
42 #include <openvdb/Platform.h>
43 #include "LevelSetTracker.h"
44 #include "Interpolation.h" // for BoxSampler, etc.
45 #include <openvdb/math/FiniteDifference.h>
46 #include <boost/math/constants/constants.hpp>
47 
48 namespace openvdb {
50 namespace OPENVDB_VERSION_NAME {
51 namespace tools {
52 
59 
64 
67 template <typename VelGridT, typename Interpolator = BoxSampler>
69 {
70 public:
71  typedef typename VelGridT::ValueType VectorType;
72  typedef typename VectorType::ValueType ScalarType;
73 
74  DiscreteField(const VelGridT &vel): mAccessor(vel.tree()), mTransform(&vel.transform()) {}
75 
79  const math::Transform& transform() const { return *mTransform; }
80 
82  inline VectorType operator() (const Vec3d& xyz, ScalarType) const
83  {
84  VectorType result = zeroVal<VectorType>();
85  Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz), result);
86  return result;
87  }
88 
90  inline VectorType operator() (const Coord& ijk, ScalarType) const
91  {
92  return mAccessor.getValue(ijk);
93  }
94 
95 private:
96  const typename VelGridT::ConstAccessor mAccessor;//Not thread-safe
97  const math::Transform* mTransform;
98 
99 }; // end of DiscreteField
100 
107 template <typename ScalarT = float>
109 {
110 public:
111  typedef ScalarT ScalarType;
113 
115 
120 
123  inline VectorType operator() (const Vec3d& xyz, ScalarType time) const;
124 
126  inline VectorType operator() (const Coord& ijk, ScalarType time) const
127  {
128  return (*this)(ijk.asVec3d(), time);
129  }
130 }; // end of EnrightField
131 
171 
172 template<typename GridT,
173  typename FieldT = EnrightField<typename GridT::ValueType>,
174  typename InterruptT = util::NullInterrupter>
176 {
177 public:
178  typedef GridT GridType;
180  typedef typename TrackerT::RangeType RangeType;
181  typedef typename TrackerT::LeafType LeafType;
183  typedef typename TrackerT::ValueType ScalarType;
184  typedef typename FieldT::VectorType VectorType;
185 
187  LevelSetAdvection(GridT& grid, const FieldT& field, InterruptT* interrupt = NULL):
188  mTracker(grid, interrupt), mField(field),
189  mSpatialScheme(math::HJWENO5_BIAS),
190  mTemporalScheme(math::TVD_RK2) {}
191 
192  virtual ~LevelSetAdvection() {};
193 
195  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
197  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
198 
200  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
202  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
203 
205  math::BiasedGradientScheme getTrackerSpatialScheme() const { return mTracker.getSpatialScheme(); }
207  void setTrackerSpatialScheme(math::BiasedGradientScheme scheme) { mTracker.setSpatialScheme(scheme); }
208 
210  math::TemporalIntegrationScheme getTrackerTemporalScheme() const { return mTracker.getTemporalScheme(); }
212  void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme) { mTracker.setTemporalScheme(scheme); }
213 
216  int getNormCount() const { return mTracker.getNormCount(); }
219  void setNormCount(int n) { mTracker.setNormCount(n); }
220 
222  int getGrainSize() const { return mTracker.getGrainSize(); }
225  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
226 
231  size_t advect(ScalarType time0, ScalarType time1);
232 
233 private:
234 
235  // This templated private class implements all the level set magic.
236  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
237  math::TemporalIntegrationScheme TemporalScheme>
238  class LevelSetAdvect
239  {
240  public:
242  LevelSetAdvect(LevelSetAdvection& parent);
244  LevelSetAdvect(const LevelSetAdvect& other);
246  LevelSetAdvect(LevelSetAdvect& other, tbb::split);
248  virtual ~LevelSetAdvect() {if (mIsMaster) this->clearField();};
251  size_t advect(ScalarType time0, ScalarType time1);
253  void operator()(const RangeType& r) const
254  {
255  if (mTask) mTask(const_cast<LevelSetAdvect*>(this), r);
256  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
257  }
259  void operator()(const RangeType& r)
260  {
261  if (mTask) mTask(this, r);
262  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
263  }
265  void join(const LevelSetAdvect& other) { mMaxAbsV = math::Max(mMaxAbsV, other.mMaxAbsV); }
266  private:
267  typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
268  LevelSetAdvection& mParent;
269  VectorType** mVec;
270  const ScalarType mMinAbsV;
271  ScalarType mMaxAbsV;
272  const MapT* mMap;
273  FuncType mTask;
274  const bool mIsMaster;
276  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
277  // method calling tbb
278  void cook(ThreadingMode mode, size_t swapBuffer = 0);
280  typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
281  void clearField();
282  void sampleXformedField(const RangeType& r, ScalarType time0, ScalarType time1);
283  void sampleAlignedField(const RangeType& r, ScalarType time0, ScalarType time1);
284  // Forward Euler advection steps: Phi(result) = Phi(0) - dt * V.Grad(0);
285  void euler1(const RangeType& r, ScalarType dt, Index resultBuffer);
286  // Convex combination of Phi and a forward Euler advection steps:
287  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * V.Grad(0));
288  void euler2(const RangeType& r, ScalarType dt, ScalarType alpha, Index phiBuffer, Index resultBuffer);
289  }; // end of private LevelSetAdvect class
290 
291  template<math::BiasedGradientScheme SpatialScheme>
292  size_t advect1(ScalarType time0, ScalarType time1);
293 
294  template<math::BiasedGradientScheme SpatialScheme,
295  math::TemporalIntegrationScheme TemporalScheme>
296  size_t advect2(ScalarType time0, ScalarType time1);
297 
298  template<math::BiasedGradientScheme SpatialScheme,
299  math::TemporalIntegrationScheme TemporalScheme,
300  typename MapType>
301  size_t advect3(ScalarType time0, ScalarType time1);
302 
303  TrackerT mTracker;
304  //each thread needs a deep copy of the field since it might contain a ValueAccessor
305  const FieldT mField;
306  math::BiasedGradientScheme mSpatialScheme;
307  math::TemporalIntegrationScheme mTemporalScheme;
308 
309  // disallow copy by assignment
310  void operator=(const LevelSetAdvection& other) {}
311 
312 };//end of LevelSetAdvection
313 
314 template<typename GridT, typename FieldT, typename InterruptT>
315 inline size_t
316 LevelSetAdvection<GridT, FieldT, InterruptT>::advect(ScalarType time0, ScalarType time1)
317 {
318  switch (mSpatialScheme) {
319  case math::FIRST_BIAS:
320  return this->advect1<math::FIRST_BIAS >(time0, time1);
321  case math::SECOND_BIAS:
322  return this->advect1<math::SECOND_BIAS >(time0, time1);
323  case math::THIRD_BIAS:
324  return this->advect1<math::THIRD_BIAS >(time0, time1);
325  case math::WENO5_BIAS:
326  return this->advect1<math::WENO5_BIAS >(time0, time1);
327  case math::HJWENO5_BIAS:
328  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
329  default:
330  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
331  }
332  return 0;
333 }
334 
335 template<typename GridT, typename FieldT, typename InterruptT>
336 template<math::BiasedGradientScheme SpatialScheme>
337 inline size_t
338 LevelSetAdvection<GridT, FieldT, InterruptT>::advect1(ScalarType time0, ScalarType time1)
339 {
340  switch (mTemporalScheme) {
341  case math::TVD_RK1:
342  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
343  case math::TVD_RK2:
344  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
345  case math::TVD_RK3:
346  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
347  default:
348  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
349  }
350  return 0;
351 }
352 
353 template<typename GridT, typename FieldT, typename InterruptT>
354 template<math::BiasedGradientScheme SpatialScheme,
355  math::TemporalIntegrationScheme TemporalScheme>
356 inline size_t
357 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
358 {
359  const math::Transform& trans = mTracker.grid().transform();
360  if (trans.mapType() == math::UniformScaleMap::mapType()) {
361  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
362  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
363  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
364  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
365  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
366  } else if (trans.mapType() == math::TranslationMap::mapType()) {
367  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
368  } else {
369  OPENVDB_THROW(ValueError, "MapType not supported!");
370  }
371  return 0;
372 }
373 
374 template<typename GridT, typename FieldT, typename InterruptT>
375 template<math::BiasedGradientScheme SpatialScheme,
376  math::TemporalIntegrationScheme TemporalScheme,
377  typename MapT>
378 inline size_t
379 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
380 {
381  LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*this);
382  return tmp.advect(time0, time1);
383 }
384 
386 
387 template <typename ScalarT>
388 inline math::Vec3<ScalarT>
389 EnrightField<ScalarT>::operator() (const Vec3d& xyz, ScalarType time) const
390 {
391  const ScalarT pi = boost::math::constants::pi<ScalarT>();
392  const ScalarT phase = pi / ScalarT(3.0);
393  const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
394  const ScalarT tr = cos(ScalarT(time) * phase);
395  const ScalarT a = sin(ScalarT(2.0)*Py);
396  const ScalarT b = -sin(ScalarT(2.0)*Px);
397  const ScalarT c = sin(ScalarT(2.0)*Pz);
398  return math::Vec3<ScalarT>(
399  tr * ( ScalarT(2) * math::Pow2(sin(Px)) * a * c ),
400  tr * ( b * math::Pow2(sin(Py)) * c ),
401  tr * ( b * a * math::Pow2(sin(Pz)) ));
402 }
403 
404 
406 
407 
408 template<typename GridT, typename FieldT, typename InterruptT>
409 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
410  math::TemporalIntegrationScheme TemporalScheme>
411 inline
415  mParent(parent),
416  mVec(NULL),
417  mMinAbsV(ScalarType(1e-6)),
418  mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()),
419  mTask(0),
420  mIsMaster(true)
421 {
422 }
423 
424 template<typename GridT, typename FieldT, typename InterruptT>
425 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
426  math::TemporalIntegrationScheme TemporalScheme>
427 inline
428 LevelSetAdvection<GridT, FieldT, InterruptT>::
429 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
430 LevelSetAdvect(const LevelSetAdvect& other):
431  mParent(other.mParent),
432  mVec(other.mVec),
433  mMinAbsV(other.mMinAbsV),
434  mMaxAbsV(other.mMaxAbsV),
435  mMap(other.mMap),
436  mTask(other.mTask),
437  mIsMaster(false)
438 {
439 }
440 
441 template<typename GridT, typename FieldT, typename InterruptT>
442 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
443  math::TemporalIntegrationScheme TemporalScheme>
444 inline
445 LevelSetAdvection<GridT, FieldT, InterruptT>::
446 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
447 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
448  mParent(other.mParent),
449  mVec(other.mVec),
450  mMinAbsV(other.mMinAbsV),
451  mMaxAbsV(other.mMaxAbsV),
452  mMap(other.mMap),
453  mTask(other.mTask),
454  mIsMaster(false)
455 {
456 }
457 
458 template<typename GridT, typename FieldT, typename InterruptT>
459 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
460  math::TemporalIntegrationScheme TemporalScheme>
461 inline size_t
464 advect(ScalarType time0, ScalarType time1)
465 {
466  size_t countCFL = 0;
467  if ( math::isZero(time0 - time1) ) return countCFL;
468  const bool isForward = time0 < time1;
469  while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
471  mParent.mTracker.leafs().rebuildAuxBuffers(TemporalScheme == math::TVD_RK3 ? 2 : 1);
472 
473  const ScalarType dt = this->sampleField(time0, time1);
474  if ( math::isZero(dt) ) break;//V is essentially zero so terminate
475 
476  OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time
477  switch(TemporalScheme) {
478  case math::TVD_RK1:
479  // Perform one explicit Euler step: t1 = t0 + dt
480  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
481  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
482  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
483  this->cook(PARALLEL_FOR, 1);
484  break;
485  case math::TVD_RK2:
486  // Perform one explicit Euler step: t1 = t0 + dt
487  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
488  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
489  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
490  this->cook(PARALLEL_FOR, 1);
491 
492  // Convex combine explict Euler step: t2 = t0 + dt
493  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * V.Grad_t1(0))
494  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), /*phi=*/1, /*result=*/1);
495  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
496  this->cook(PARALLEL_FOR, 1);
497  break;
498  case math::TVD_RK3:
499  // Perform one explicit Euler step: t1 = t0 + dt
500  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
501  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
502  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
503  this->cook(PARALLEL_FOR, 1);
504 
505  // Convex combine explict Euler step: t2 = t0 + dt/2
506  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * V.Grad_t1(0))
507  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), /*phi=*/1, /*result=*/2);
508  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
509  this->cook(PARALLEL_FOR, 2);
510 
511  // Convex combine explict Euler step: t3 = t0 + dt
512  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * V.Grad_t2(0)
513  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), /*phi=*/1, /*result=*/2);
514  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
515  this->cook(PARALLEL_FOR, 2);
516  break;
517  default:
518  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
519  }//end of compile-time resolved switch
521 
522  time0 += isForward ? dt : -dt;
523  ++countCFL;
524  mParent.mTracker.leafs().removeAuxBuffers();
525  this->clearField();
527  mParent.mTracker.track();
528  }//end wile-loop over time
529  return countCFL;//number of CLF propagation steps
530 }
531 
532 template<typename GridT, typename FieldT, typename InterruptT>
533 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
534  math::TemporalIntegrationScheme TemporalScheme>
535 inline typename GridT::ValueType
536 LevelSetAdvection<GridT, FieldT, InterruptT>::
537 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
538 sampleField(ScalarType time0, ScalarType time1)
539 {
540  mMaxAbsV = mMinAbsV;
541  const size_t leafCount = mParent.mTracker.leafs().leafCount();
542  if (leafCount==0) return ScalarType(0.0);
543  mVec = new VectorType*[leafCount];
544  if (mParent.mField.transform() == mParent.mTracker.grid().transform()) {
545  mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0, time1);
546  } else {
547  mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0, time1);
548  }
549  this->cook(PARALLEL_REDUCE);
550  if (math::isExactlyEqual(mMinAbsV, mMaxAbsV)) return ScalarType(0.0);//V is essentially zero
551 #ifndef _MSC_VER // Visual C++ doesn't guarantee thread-safe initialization of local statics
552  static
553 #endif
554  const ScalarType CFL = (TemporalScheme == math::TVD_RK1 ? ScalarType(0.3) :
555  TemporalScheme == math::TVD_RK2 ? ScalarType(0.9) :
556  ScalarType(1.0))/math::Sqrt(ScalarType(3.0));
557  const ScalarType dt = math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
558  return math::Min(dt, ScalarType(CFL*dx/math::Sqrt(mMaxAbsV)));
559 }
560 
561 template<typename GridT, typename FieldT, typename InterruptT>
562 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
563  math::TemporalIntegrationScheme TemporalScheme>
564 inline void
565 LevelSetAdvection<GridT, FieldT, InterruptT>::
566 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
567 sampleXformedField(const RangeType& range, ScalarType time0, ScalarType time1)
568 {
569  const bool isForward = time0 < time1;
570  typedef typename LeafType::ValueOnCIter VoxelIterT;
571  const MapT& map = *mMap;
572  mParent.mTracker.checkInterrupter();
573  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
574  const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
575  VectorType* vec = new VectorType[leaf.onVoxelCount()];
576  int m = 0;
577  for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
578  const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
579  mMaxAbsV = math::Max(mMaxAbsV, ScalarType(math::Pow2(V[0])+math::Pow2(V[1])+math::Pow2(V[2])));
580  vec[m] = isForward ? V : -V;
581  }
582  mVec[n] = vec;
583  }
584 }
585 
586 template<typename GridT, typename FieldT, typename InterruptT>
587 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
588  math::TemporalIntegrationScheme TemporalScheme>
589 inline void
590 LevelSetAdvection<GridT, FieldT, InterruptT>::
591 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
592 sampleAlignedField(const RangeType& range, ScalarType time0, ScalarType time1)
593 {
594  const bool isForward = time0 < time1;
595  typedef typename LeafType::ValueOnCIter VoxelIterT;
596  mParent.mTracker.checkInterrupter();
597  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
598  const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
599  VectorType* vec = new VectorType[leaf.onVoxelCount()];
600  int m = 0;
601  for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
602  const VectorType V = mParent.mField(iter.getCoord(), time0);
603  mMaxAbsV = math::Max(mMaxAbsV, ScalarType(math::Pow2(V[0])+math::Pow2(V[1])+math::Pow2(V[2])));
604  vec[m] = isForward ? V : -V;
605  }
606  mVec[n] = vec;
607  }
608 }
609 
610 template<typename GridT, typename FieldT, typename InterruptT>
611 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
612  math::TemporalIntegrationScheme TemporalScheme>
613 inline void
614 LevelSetAdvection<GridT, FieldT, InterruptT>::
615 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
616 clearField()
617 {
618  if (mVec == NULL) return;
619  for (size_t n=0, e=mParent.mTracker.leafs().leafCount(); n<e; ++n) delete [] mVec[n];
620  delete [] mVec;
621  mVec = NULL;
622 }
623 
624 template<typename GridT, typename FieldT, typename InterruptT>
625 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
626  math::TemporalIntegrationScheme TemporalScheme>
627 inline void
628 LevelSetAdvection<GridT, FieldT, InterruptT>::
629 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
630 cook(ThreadingMode mode, size_t swapBuffer)
631 {
632  mParent.mTracker.startInterrupter("Advecting level set");
633 
634  if (mParent.mTracker.getGrainSize()==0) {
635  (*this)(mParent.mTracker.leafs().getRange());
636  } else if (mode == PARALLEL_FOR) {
637  tbb::parallel_for(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *this);
638  } else if (mode == PARALLEL_REDUCE) {
639  tbb::parallel_reduce(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *this);
640  } else {
641  throw std::runtime_error("Undefined threading mode");
642  }
643 
644  mParent.mTracker.leafs().swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
645 
646  mParent.mTracker.endInterrupter();
647 }
648 
649 // Forward Euler advection steps:
650 // Phi(result) = Phi(0) - dt * V.Grad(0);
651 template<typename GridT, typename FieldT, typename InterruptT>
652 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
653  math::TemporalIntegrationScheme TemporalScheme>
654 inline void
655 LevelSetAdvection<GridT, FieldT, InterruptT>::
656 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
657 euler1(const RangeType& range, ScalarType dt, Index resultBuffer)
658 {
659  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
660  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
661  typedef typename LeafType::ValueOnCIter VoxelIterT;
662  mParent.mTracker.checkInterrupter();
663  const MapT& map = *mMap;
664  typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
665  Stencil stencil(mParent.mTracker.grid());
666  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
667  BufferType& result = leafs.getBuffer(n, resultBuffer);
668  const VectorType* vec = mVec[n];
669  int m=0;
670  for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
671  stencil.moveTo(iter);
672  const VectorType V = vec[m], G = math::GradientBiased<MapT, SpatialScheme>::result(map, stencil, V);
673  result.setValue(iter.pos(), *iter - dt * V.dot(G));
674  }
675  }
676 }
677 
678 // Convex combination of Phi and a forward Euler advection steps:
679 // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * V.Grad(0));
680 template<typename GridT, typename FieldT, typename InterruptT>
681 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
682  math::TemporalIntegrationScheme TemporalScheme>
683 inline void
684 LevelSetAdvection<GridT, FieldT, InterruptT>::
685 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
686 euler2(const RangeType& range, ScalarType dt, ScalarType alpha, Index phiBuffer, Index resultBuffer)
687 {
688  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
689  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
690  typedef typename LeafType::ValueOnCIter VoxelIterT;
691  mParent.mTracker.checkInterrupter();
692  const MapT& map = *mMap;
693  typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
694  const ScalarType beta = ScalarType(1.0) - alpha;
695  Stencil stencil(mParent.mTracker.grid());
696  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
697  const BufferType& phi = leafs.getBuffer(n, phiBuffer);
698  BufferType& result = leafs.getBuffer(n, resultBuffer);
699  const VectorType* vec = mVec[n];
700  int m=0;
701  for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
702  stencil.moveTo(iter);
703  const VectorType V = vec[m], G = math::GradientBiased<MapT, SpatialScheme>::result(map, stencil, V);
704  result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
705  }
706  }
707 }
708 
709 } // namespace tools
710 } // namespace OPENVDB_VERSION_NAME
711 } // namespace openvdb
712 
713 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
714 
715 // Copyright (c) 2012-2014 DreamWorks Animation LLC
716 // All rights reserved. This software is distributed under the
717 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
ScalarT ScalarType
Definition: LevelSetAdvect.h:111
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:265
Thin wrapper class for a velocity grid.
Definition: LevelSetAdvect.h:68
math::Vec3< ScalarT > VectorType
Definition: LevelSetAdvect.h:112
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite difference scheme.
Definition: LevelSetAdvect.h:197
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:284
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
Name mapType() const
Return the transformation map's type-name.
Definition: Transform.h:93
Definition: Mat.h:146
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
TrackerT::BufferType BufferType
Definition: LevelSetAdvect.h:182
TrackerT::ValueType ScalarType
Definition: LevelSetAdvect.h:183
math::BiasedGradientScheme getTrackerSpatialScheme() const
Definition: LevelSetAdvect.h:205
VelGridT::ValueType VectorType
Definition: LevelSetAdvect.h:71
static math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk, const Vec3< typename Accessor::ValueType > &V)
Definition: Operators.h:831
DiscreteField(const VelGridT &vel)
Definition: LevelSetAdvect.h:74
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
LevelSetTracker< GridT, InterruptT > TrackerT
Definition: LevelSetAdvect.h:179
Type Pow2(Type x)
Return .
Definition: Math.h:498
TreeType::LeafNodeType LeafType
Definition: LevelSetTracker.h:73
TrackerT::LeafType LeafType
Definition: LevelSetAdvect.h:181
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:606
const math::Transform & transform() const
Definition: LevelSetAdvect.h:79
virtual ~LevelSetAdvection()
Definition: LevelSetAdvect.h:192
math::Transform transform() const
Definition: LevelSetAdvect.h:119
Definition: FiniteDifference.h:196
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the spatial finite difference scheme.
Definition: LevelSetAdvect.h:212
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the spatial finite difference scheme.
Definition: LevelSetAdvect.h:202
size_t advect(ScalarType time0, ScalarType time1)
Definition: LevelSetAdvect.h:316
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:693
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Index32 Index
Definition: Types.h:59
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite difference scheme.
Definition: LevelSetAdvect.h:207
LevelSetAdvection(GridT &grid, const FieldT &field, InterruptT *interrupt=NULL)
Main constructor.
Definition: LevelSetAdvect.h:187
int getNormCount() const
Definition: LevelSetAdvect.h:216
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Definition: LevelSetAdvect.h:210
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
Definition: LevelSetAdvect.h:225
Definition: Exceptions.h:39
Definition: FiniteDifference.h:195
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:121
Vec3< double > Vec3d
Definition: Vec3.h:629
math::TemporalIntegrationScheme getTemporalScheme() const
Definition: LevelSetAdvect.h:200
LeafManagerType::RangeType RangeType
Definition: LevelSetTracker.h:76
math::BiasedGradientScheme getSpatialScheme() const
Definition: LevelSetAdvect.h:195
Definition: FiniteDifference.h:264
VectorType::ValueType ScalarType
Definition: LevelSetAdvect.h:72
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
Definition: LevelSetAdvect.h:219
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:545
Definition: Exceptions.h:88
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
FieldT::VectorType VectorType
Definition: LevelSetAdvect.h:184
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:122
GridT GridType
Definition: LevelSetAdvect.h:178
Hyperbolic advection of narrow-band level sets in an external velocity field.
Definition: LevelSetAdvect.h:175
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:391
int getGrainSize() const
Definition: LevelSetAdvect.h:222
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:68
TrackerT::RangeType RangeType
Definition: LevelSetAdvect.h:180
Definition: FiniteDifference.h:198
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
LeafManagerType::BufferType BufferType
Definition: LevelSetTracker.h:78
Definition: FiniteDifference.h:194
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:314
EnrightField()
Definition: LevelSetAdvect.h:114
Definition: FiniteDifference.h:263
TreeType::ValueType ValueType
Definition: LevelSetTracker.h:74
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Analytical, divergence-free and periodic vecloity field.
Definition: LevelSetAdvect.h:108