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