OpenVDB  1.2.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 
101 template <typename ScalarT = float>
103 {
104 public:
105  typedef ScalarT ScalarType;
107 
109 
114 
117  inline VectorType operator() (const Vec3d& xyz, ScalarType time) const;
118 
120  inline VectorType operator() (const Coord& ijk, ScalarType time) const
121  {
122  return (*this)(ijk.asVec3d(), time);
123  }
124 }; // end of EnrightField
125 
165 
166 template<typename GridT,
167  typename FieldT = EnrightField<typename GridT::ValueType>,
168  typename InterruptT = util::NullInterrupter>
170 {
171 public:
172  typedef GridT GridType;
174  typedef typename TrackerT::RangeType RangeType;
175  typedef typename TrackerT::LeafType LeafType;
177  typedef typename TrackerT::ValueType ScalarType;
178  typedef typename FieldT::VectorType VectorType;
179 
181  LevelSetAdvection(GridT& grid, const FieldT& field, InterruptT* interrupt = NULL):
182  mTracker(grid, interrupt), mField(field),
183  mSpatialScheme(math::HJWENO5_BIAS),
184  mTemporalScheme(math::TVD_RK2) {}
185 
186  virtual ~LevelSetAdvection() {};
187 
189  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
191  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
192 
194  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
196  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
197 
199  math::BiasedGradientScheme getTrackerSpatialScheme() const { return mTracker.getSpatialScheme(); }
201  void setTrackerSpatialScheme(math::BiasedGradientScheme scheme) { mTracker.setSpatialScheme(scheme); }
202 
204  math::TemporalIntegrationScheme getTrackerTemporalScheme() const { return mTracker.getTemporalScheme(); }
206  void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme) { mTracker.setTemporalScheme(scheme); }
207 
210  int getNormCount() const { return mTracker.getNormCount(); }
213  void setNormCount(int n) { mTracker.setNormCount(n); }
214 
216  int getGrainSize() const { return mTracker.getGrainSize(); }
219  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
220 
225  size_t advect(ScalarType time0, ScalarType time1);
226 
227 private:
228 
229  // This templated private class implements all the level set magic.
230  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
231  math::TemporalIntegrationScheme TemporalScheme>
232  class LevelSetAdvect
233  {
234  public:
236  LevelSetAdvect(LevelSetAdvection& parent);
238  LevelSetAdvect(const LevelSetAdvect& other);
240  LevelSetAdvect(LevelSetAdvect& other, tbb::split);
242  virtual ~LevelSetAdvect() {if (mIsMaster) this->clearField();};
245  size_t advect(ScalarType time0, ScalarType time1);
247  void operator()(const RangeType& r) const
248  {
249  if (mTask) mTask(const_cast<LevelSetAdvect*>(this), r);
250  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
251  }
253  void operator()(const RangeType& r)
254  {
255  if (mTask) mTask(this, r);
256  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
257  }
259  void join(const LevelSetAdvect& other) { mMaxAbsV = math::Max(mMaxAbsV, other.mMaxAbsV); }
260  private:
261  typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
262  LevelSetAdvection& mParent;
263  VectorType** mVec;
264  const ScalarType mMinAbsV;
265  ScalarType mMaxAbsV;
266  const MapT* mMap;
267  FuncType mTask;
268  const bool mIsMaster;
270  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
271  // method calling tbb
272  void cook(ThreadingMode mode, size_t swapBuffer = 0);
274  typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
275  void clearField();
276  void sampleXformedField(const RangeType& r, ScalarType time0, ScalarType time1);
277  void sampleAlignedField(const RangeType& r, ScalarType time0, ScalarType time1);
278  // Forward Euler advection steps: Phi(result) = Phi(0) - dt * V.Grad(0);
279  void euler1(const RangeType& r, ScalarType dt, Index resultBuffer);
280  // Convex combination of Phi and a forward Euler advection steps:
281  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * V.Grad(0));
282  void euler2(const RangeType& r, ScalarType dt, ScalarType alpha, Index phiBuffer, Index resultBuffer);
283  }; // end of private LevelSetAdvect class
284 
285  template<math::BiasedGradientScheme SpatialScheme>
286  size_t advect1(ScalarType time0, ScalarType time1);
287 
288  template<math::BiasedGradientScheme SpatialScheme,
289  math::TemporalIntegrationScheme TemporalScheme>
290  size_t advect2(ScalarType time0, ScalarType time1);
291 
292  template<math::BiasedGradientScheme SpatialScheme,
293  math::TemporalIntegrationScheme TemporalScheme,
294  typename MapType>
295  size_t advect3(ScalarType time0, ScalarType time1);
296 
297  TrackerT mTracker;
298  //each thread needs a deep copy of the field since it might contain a ValueAccessor
299  const FieldT mField;
300  math::BiasedGradientScheme mSpatialScheme;
301  math::TemporalIntegrationScheme mTemporalScheme;
302 
303  // disallow copy by assignment
304  void operator=(const LevelSetAdvection& other) {}
305 
306 };//end of LevelSetAdvection
307 
308 template<typename GridT, typename FieldT, typename InterruptT>
309 inline size_t
311 {
312  switch (mSpatialScheme) {
313  case math::FIRST_BIAS:
314  return this->advect1<math::FIRST_BIAS >(time0, time1);
315  case math::SECOND_BIAS:
316  return this->advect1<math::SECOND_BIAS >(time0, time1);
317  case math::THIRD_BIAS:
318  return this->advect1<math::THIRD_BIAS >(time0, time1);
319  case math::WENO5_BIAS:
320  return this->advect1<math::WENO5_BIAS >(time0, time1);
321  case math::HJWENO5_BIAS:
322  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
323  default:
324  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
325  }
326  return 0;
327 }
328 
329 template<typename GridT, typename FieldT, typename InterruptT>
330 template<math::BiasedGradientScheme SpatialScheme>
331 inline size_t
332 LevelSetAdvection<GridT, FieldT, InterruptT>::advect1(ScalarType time0, ScalarType time1)
333 {
334  switch (mTemporalScheme) {
335  case math::TVD_RK1:
336  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
337  case math::TVD_RK2:
338  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
339  case math::TVD_RK3:
340  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
341  default:
342  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
343  }
344  return 0;
345 }
346 
347 template<typename GridT, typename FieldT, typename InterruptT>
348 template<math::BiasedGradientScheme SpatialScheme,
349  math::TemporalIntegrationScheme TemporalScheme>
350 inline size_t
351 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
352 {
353  const math::Transform& trans = mTracker.grid().transform();
354  if (trans.mapType() == math::UniformScaleMap::mapType()) {
355  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
356  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
357  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
358  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
359  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
360  } else if (trans.mapType() == math::TranslationMap::mapType()) {
361  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
362  } else {
363  OPENVDB_THROW(ValueError, "MapType not supported!");
364  }
365  return 0;
366 }
367 
368 template<typename GridT, typename FieldT, typename InterruptT>
369 template<math::BiasedGradientScheme SpatialScheme,
370  math::TemporalIntegrationScheme TemporalScheme,
371  typename MapT>
372 inline size_t
373 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
374 {
375  LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*this);
376  return tmp.advect(time0, time1);
377 }
378 
380 
381 template <typename ScalarT>
382 inline math::Vec3<ScalarT>
384 {
385  static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
386  const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
387  const ScalarT tr = cos(ScalarT(time) * phase);
388  const ScalarT a = sin(ScalarT(2.0)*Py);
389  const ScalarT b = -sin(ScalarT(2.0)*Px);
390  const ScalarT c = sin(ScalarT(2.0)*Pz);
391  return math::Vec3<ScalarT>(
392  tr * ( ScalarT(2) * math::Pow2(sin(Px)) * a * c ),
393  tr * ( b * math::Pow2(sin(Py)) * c ),
394  tr * ( b * a * math::Pow2(sin(Pz)) ));
395 }
396 
397 
399 
400 
401 template<typename GridT, typename FieldT, typename InterruptT>
402 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
403  math::TemporalIntegrationScheme TemporalScheme>
404 inline
408  mParent(parent),
409  mVec(NULL),
410  mMinAbsV(1e-6),
411  mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()),
412  mTask(0),
413  mIsMaster(true)
414 {
415 }
416 
417 template<typename GridT, typename FieldT, typename InterruptT>
418 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
419  math::TemporalIntegrationScheme TemporalScheme>
420 inline
421 LevelSetAdvection<GridT, FieldT, InterruptT>::
422 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
423 LevelSetAdvect(const LevelSetAdvect& other):
424  mParent(other.mParent),
425  mVec(other.mVec),
426  mMinAbsV(other.mMinAbsV),
427  mMaxAbsV(other.mMaxAbsV),
428  mMap(other.mMap),
429  mTask(other.mTask),
430  mIsMaster(false)
431 {
432 }
433 
434 template<typename GridT, typename FieldT, typename InterruptT>
435 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
436  math::TemporalIntegrationScheme TemporalScheme>
437 inline
438 LevelSetAdvection<GridT, FieldT, InterruptT>::
439 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
440 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
441  mParent(other.mParent),
442  mVec(other.mVec),
443  mMinAbsV(other.mMinAbsV),
444  mMaxAbsV(other.mMaxAbsV),
445  mMap(other.mMap),
446  mTask(other.mTask),
447  mIsMaster(false)
448 {
449 }
450 
451 template<typename GridT, typename FieldT, typename InterruptT>
452 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
453  math::TemporalIntegrationScheme TemporalScheme>
454 inline size_t
457 advect(ScalarType time0, ScalarType time1)
458 {
459  size_t countCFL = 0;
460  if ( math::isZero(time0 - time1) ) return countCFL;
461  const bool isForward = time0 < time1;
462  while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
464  mParent.mTracker.leafs().rebuildAuxBuffers(TemporalScheme == math::TVD_RK3 ? 2 : 1);
465 
466  const ScalarType dt = this->sampleField(time0, time1);
467  if ( math::isZero(dt) ) break;//V is essentially zero so terminate
468 
469  switch(TemporalScheme) {//switch is resolved at compile-time
470  case math::TVD_RK1:
471  // Perform one explicit Euler step: t1 = t0 + dt
472  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
473  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
474  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
475  this->cook(PARALLEL_FOR, 1);
476  break;
477  case math::TVD_RK2:
478  // Perform one explicit Euler step: t1 = t0 + dt
479  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
480  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
481  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
482  this->cook(PARALLEL_FOR, 1);
483 
484  // Convex combine explict Euler step: t2 = t0 + dt
485  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * V.Grad_t1(0))
486  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), /*phi=*/1, /*result=*/1);
487  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
488  this->cook(PARALLEL_FOR, 1);
489  break;
490  case math::TVD_RK3:
491  // Perform one explicit Euler step: t1 = t0 + dt
492  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
493  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
494  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
495  this->cook(PARALLEL_FOR, 1);
496 
497  // Convex combine explict Euler step: t2 = t0 + dt/2
498  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * V.Grad_t1(0))
499  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), /*phi=*/1, /*result=*/2);
500  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
501  this->cook(PARALLEL_FOR, 2);
502 
503  // Convex combine explict Euler step: t3 = t0 + dt
504  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * V.Grad_t2(0)
505  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), /*phi=*/1, /*result=*/2);
506  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
507  this->cook(PARALLEL_FOR, 2);
508  break;
509  default:
510  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
511  }
512  time0 += isForward ? dt : -dt;
513  ++countCFL;
514  mParent.mTracker.leafs().removeAuxBuffers();
515  this->clearField();
517  mParent.mTracker.track();
518  }//end wile-loop over time
519  return countCFL;//number of CLF propagation steps
520 }
521 
522 template<typename GridT, typename FieldT, typename InterruptT>
523 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
524  math::TemporalIntegrationScheme TemporalScheme>
525 inline typename GridT::ValueType
526 LevelSetAdvection<GridT, FieldT, InterruptT>::
527 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
528 sampleField(ScalarType time0, ScalarType time1)
529 {
530  mMaxAbsV = mMinAbsV;
531  const size_t leafCount = mParent.mTracker.leafs().leafCount();
532  if (leafCount==0) return ScalarType(0.0);
533  mVec = new VectorType*[leafCount];
534  if (mParent.mField.transform() == mParent.mTracker.grid().transform()) {
535  mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0, time1);
536  } else {
537  mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0, time1);
538  }
539  this->cook(PARALLEL_REDUCE);
540  if (math::isExactlyEqual(mMinAbsV, mMaxAbsV)) return ScalarType(0.0);//V is essentially zero
541  static const ScalarType CFL = (TemporalScheme == math::TVD_RK1 ? ScalarType(0.3) :
542  TemporalScheme == math::TVD_RK2 ? ScalarType(0.9) :
543  ScalarType(1.0))/math::Sqrt(ScalarType(3.0));
544  const ScalarType dt = math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
545  return math::Min(dt, ScalarType(CFL*dx/math::Sqrt(mMaxAbsV)));
546 }
547 
548 template<typename GridT, typename FieldT, typename InterruptT>
549 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
550  math::TemporalIntegrationScheme TemporalScheme>
551 inline void
552 LevelSetAdvection<GridT, FieldT, InterruptT>::
553 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
554 sampleXformedField(const RangeType& range, ScalarType time0, ScalarType time1)
555 {
556  const bool isForward = time0 < time1;
557  typedef typename LeafType::ValueOnCIter VoxelIterT;
558  const MapT& map = *mMap;
559  mParent.mTracker.checkInterrupter();
560  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
561  const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
562  VectorType* vec = new VectorType[leaf.onVoxelCount()];
563  int m = 0;
564  for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
565  const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
566  mMaxAbsV = math::Max(mMaxAbsV, ScalarType(math::Pow2(V[0])+math::Pow2(V[1])+math::Pow2(V[2])));
567  vec[m] = isForward ? V : -V;
568  }
569  mVec[n] = vec;
570  }
571 }
572 
573 template<typename GridT, typename FieldT, typename InterruptT>
574 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
575  math::TemporalIntegrationScheme TemporalScheme>
576 inline void
577 LevelSetAdvection<GridT, FieldT, InterruptT>::
578 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
579 sampleAlignedField(const RangeType& range, ScalarType time0, ScalarType time1)
580 {
581  const bool isForward = time0 < time1;
582  typedef typename LeafType::ValueOnCIter VoxelIterT;
583  mParent.mTracker.checkInterrupter();
584  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
585  const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
586  VectorType* vec = new VectorType[leaf.onVoxelCount()];
587  int m = 0;
588  for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
589  const VectorType V = mParent.mField(iter.getCoord(), time0);
590  mMaxAbsV = math::Max(mMaxAbsV, ScalarType(math::Pow2(V[0])+math::Pow2(V[1])+math::Pow2(V[2])));
591  vec[m] = isForward ? V : -V;
592  }
593  mVec[n] = vec;
594  }
595 }
596 
597 template<typename GridT, typename FieldT, typename InterruptT>
598 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
599  math::TemporalIntegrationScheme TemporalScheme>
600 inline void
601 LevelSetAdvection<GridT, FieldT, InterruptT>::
602 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
603 clearField()
604 {
605  if (mVec == NULL) return;
606  for (size_t n=0, e=mParent.mTracker.leafs().leafCount(); n<e; ++n) delete [] mVec[n];
607  delete [] mVec;
608  mVec = NULL;
609 }
610 
611 template<typename GridT, typename FieldT, typename InterruptT>
612 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
613  math::TemporalIntegrationScheme TemporalScheme>
614 inline void
615 LevelSetAdvection<GridT, FieldT, InterruptT>::
616 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
617 cook(ThreadingMode mode, size_t swapBuffer)
618 {
619  mParent.mTracker.startInterrupter("Advecting level set");
620 
621  if (mParent.mTracker.getGrainSize()==0) {
622  (*this)(mParent.mTracker.leafs().getRange());
623  } else if (mode == PARALLEL_FOR) {
624  tbb::parallel_for(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *this);
625  } else if (mode == PARALLEL_REDUCE) {
626  tbb::parallel_reduce(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *this);
627  } else {
628  throw std::runtime_error("Undefined threading mode");
629  }
630 
631  mParent.mTracker.leafs().swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
632 
633  mParent.mTracker.endInterrupter();
634 }
635 
636 // Forward Euler advection steps:
637 // Phi(result) = Phi(0) - dt * V.Grad(0);
638 template<typename GridT, typename FieldT, typename InterruptT>
639 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
640  math::TemporalIntegrationScheme TemporalScheme>
641 inline void
642 LevelSetAdvection<GridT, FieldT, InterruptT>::
643 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
644 euler1(const RangeType& range, ScalarType dt, Index resultBuffer)
645 {
646  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
647  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
648  typedef typename LeafType::ValueOnCIter VoxelIterT;
649  mParent.mTracker.checkInterrupter();
650  const MapT& map = *mMap;
651  typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
652  Stencil stencil(mParent.mTracker.grid());
653  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
654  BufferType& result = leafs.getBuffer(n, resultBuffer);
655  const VectorType* vec = mVec[n];
656  int m=0;
657  for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
658  stencil.moveTo(iter);
659  const VectorType V = vec[m], G = math::GradientBiased<MapT, SpatialScheme>::result(map, stencil, V);
660  result.setValue(iter.pos(), *iter - dt * V.dot(G));
661  }
662  }
663 }
664 
665 // Convex combination of Phi and a forward Euler advection steps:
666 // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * V.Grad(0));
667 template<typename GridT, typename FieldT, typename InterruptT>
668 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
669  math::TemporalIntegrationScheme TemporalScheme>
670 inline void
671 LevelSetAdvection<GridT, FieldT, InterruptT>::
672 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
673 euler2(const RangeType& range, ScalarType dt, ScalarType alpha, Index phiBuffer, Index resultBuffer)
674 {
675  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
676  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
677  typedef typename LeafType::ValueOnCIter VoxelIterT;
678  mParent.mTracker.checkInterrupter();
679  const MapT& map = *mMap;
680  typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
681  const ScalarType beta = ScalarType(1.0) - alpha;
682  Stencil stencil(mParent.mTracker.grid());
683  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
684  const BufferType& phi = leafs.getBuffer(n, phiBuffer);
685  BufferType& result = leafs.getBuffer(n, resultBuffer);
686  const VectorType* vec = mVec[n];
687  int m=0;
688  for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
689  stencil.moveTo(iter);
690  const VectorType V = vec[m], G = math::GradientBiased<MapT, SpatialScheme>::result(map, stencil, V);
691  result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
692  }
693  }
694 }
695 
696 } // namespace tools
697 } // namespace OPENVDB_VERSION_NAME
698 } // namespace openvdb
699 
700 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
701 
702 // Copyright (c) 2012-2013 DreamWorks Animation LLC
703 // All rights reserved. This software is distributed under the
704 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )