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