OpenVDB  0.104.0
LevelSetAdvect.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012 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 //
36 
37 #ifndef OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
38 #define OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
39 
40 #include "LevelSetTracker.h"
41 #include <openvdb/tools/GridTransformer.h>// for BoxSampler etc
42 #include <openvdb/math/FiniteDifference.h>
43 
44 namespace openvdb {
46 namespace OPENVDB_VERSION_NAME {
47 namespace tools {
48 
55 
60 
63 template <typename VelGridT, typename Interpolator = BoxSampler>
65 {
66 public:
67  typedef typename VelGridT::ConstAccessor AccessorType;
68  typedef typename VelGridT::ValueType VectorType;
69  typedef typename VectorType::ValueType ScalarType;
70 
71  DiscreteField(const VelGridT &vel):
72  mAccessor(vel.getConstAccessor()),
73  mTransform(vel.transform()) {}
74 
78  const math::Transform& transform() const { return mTransform; }
79 
81  inline VectorType operator() (const Vec3d& xyz, ScalarType) const;
82 
84  inline VectorType operator() (const Coord& ijk, ScalarType) const
85  {
86  return mAccessor.getValue(ijk);
87  }
88 
89 private:
90  const AccessorType mAccessor;
91  const math::Transform& mTransform;
92 
93 }; // end of DiscreteField
94 
97 template <typename ScalarT = float>
99 {
100 public:
101  typedef ScalarT ScalarType;
103 
105 
109  const math::Transform& transform() const { return mTransform; }
110 
113  inline VectorType operator() (const Vec3d& xyz, ScalarType time) const;
114 
116  inline VectorType operator() (const Coord& ijk, ScalarType time) const
117  {
118  return (*this)(ijk.asVec3d(), time);
119  }
120 
121 private:
122  const math::Transform mTransform;//identity transform between world and index space
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 
223  size_t advect(ScalarType time0, ScalarType time1);
224 
225 private:
226 
227  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
228  math::TemporalIntegrationScheme TemporalScheme>
229  class LevelSetAdvect
230  {
231  public:
233  LevelSetAdvect(LevelSetAdvection& parent);
235  LevelSetAdvect(const LevelSetAdvect& other);
237  LevelSetAdvect(LevelSetAdvect& other, tbb::split);
239  virtual ~LevelSetAdvect() {if (mIsMaster) this->clearField();};
242  size_t advect(ScalarType time0, ScalarType time1);
244  void operator()(const RangeType& r) const
245  {
246  if (mTask) mTask(const_cast<LevelSetAdvect*>(this), r);
247  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
248  }
250  void operator()(const RangeType& r)
251  {
252  if (mTask) mTask(this, r);
253  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
254  }
256  void join(const LevelSetAdvect& other) { mMaxAbsV = math::Max(mMaxAbsV, other.mMaxAbsV); }
257  private:
258  typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
259  LevelSetAdvection& mParent;
260  VectorType** mVec;
261  const ScalarType mMinAbsV;
262  ScalarType mMaxAbsV;
263  const MapT& mMap;
264  FuncType mTask;
265  const bool mIsMaster;
267  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
268  // method calling tbb
269  void cook(ThreadingMode mode, size_t swapBuffer = 0);
271  typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
272  void clearField();
273  void sampleXformedField(const RangeType& r, ScalarType t);
274  void sampleAlignedField(const RangeType& r, ScalarType t);
275  // Forward Euler advection steps: Phi(result) = Phi(0) - dt * V.Grad(0);
276  void euler1(const RangeType& r, ScalarType dt, Index resultBuffer);
277  // Convex combination of Phi and a forward Euler advection steps:
278  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * V.Grad(0));
279  void euler2(const RangeType& r, ScalarType dt, ScalarType alpha, Index phiBuffer, Index resultBuffer);
280  }; // end of private LevelSetAdvect class
281 
282  template<math::BiasedGradientScheme SpatialScheme>
283  size_t advect1(ScalarType time0, ScalarType time1);
284 
285  template<math::BiasedGradientScheme SpatialScheme,
286  math::TemporalIntegrationScheme TemporalScheme>
287  size_t advect2(ScalarType time0, ScalarType time1);
288 
289  template<math::BiasedGradientScheme SpatialScheme,
290  math::TemporalIntegrationScheme TemporalScheme,
291  typename MapType>
292  size_t advect3(ScalarType time0, ScalarType time1);
293 
294  TrackerT mTracker;
295  //each thread needs a deep copy of the field since it might contain a ValueAccessor
296  const FieldT mField;
297  math::BiasedGradientScheme mSpatialScheme;
298  math::TemporalIntegrationScheme mTemporalScheme;
299 
300  // disallow copy by assignment
301  void operator=(const LevelSetAdvection& other) {}
302 
303 };//end of LevelSetAdvection
304 
305 template<typename GridT, typename FieldT, typename InterruptT>
306 inline size_t
308 {
309  switch (mSpatialScheme) {
310  case math::FIRST_BIAS:
311  return this->advect1<math::FIRST_BIAS >(time0, time1);
312  case math::SECOND_BIAS:
313  return this->advect1<math::SECOND_BIAS >(time0, time1);
314  case math::THIRD_BIAS:
315  return this->advect1<math::THIRD_BIAS >(time0, time1);
316  case math::WENO5_BIAS:
317  return this->advect1<math::WENO5_BIAS >(time0, time1);
318  case math::HJWENO5_BIAS:
319  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
320  default:
321  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
322  }
323  return 0;
324 }
325 
326 template<typename GridT, typename FieldT, typename InterruptT>
327 template<math::BiasedGradientScheme SpatialScheme>
328 inline size_t
329 LevelSetAdvection<GridT, FieldT, InterruptT>::advect1(ScalarType time0, ScalarType time1)
330 {
331  switch (mTemporalScheme) {
332  case math::TVD_RK1:
333  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
334  case math::TVD_RK2:
335  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
336  case math::TVD_RK3:
337  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
338  default:
339  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
340  }
341  return 0;
342 }
343 
344 template<typename GridT, typename FieldT, typename InterruptT>
345 template<math::BiasedGradientScheme SpatialScheme,
346  math::TemporalIntegrationScheme TemporalScheme>
347 inline size_t
348 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
349 {
350  const math::Transform& trans = mTracker.grid().transform();
351  if (trans.mapType() == math::UniformScaleMap::mapType()) {
352  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
353  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
354  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
355  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
356  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
357  } else if (trans.mapType() == math::TranslationMap::mapType()) {
358  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
359  } else {
360  OPENVDB_THROW(ValueError, "MapType not supported!");
361  }
362  return 0;
363 }
364 
365 template<typename GridT, typename FieldT, typename InterruptT>
366 template<math::BiasedGradientScheme SpatialScheme,
367  math::TemporalIntegrationScheme TemporalScheme,
368  typename MapT>
369 inline size_t
370 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
371 {
372  LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*this);
373  return tmp.advect(time0, time1);
374 }
375 
376 template <typename VelGridT, typename Interpolator>
377 inline typename VelGridT::ValueType
379 {
380  const VectorType coord = mTransform.worldToIndex(xyz);
381  VectorType result;
382  Interpolator::template sample<AccessorType>(mAccessor, coord, result);
383  return result;
384 }
385 
386 template <typename ScalarT>
387 inline math::Vec3<ScalarT>
389 {
390  static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
391  const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
392  const ScalarT tr = cos(ScalarT(time) * phase);
393  const ScalarT a = sin(ScalarT(2.0)*Py);
394  const ScalarT b = -sin(ScalarT(2.0)*Px);
395  const ScalarT c = sin(ScalarT(2.0)*Pz);
396  return math::Vec3<ScalarT>(
397  tr * ( ScalarT(2) * math::Pow2(sin(Px)) * a * c ),
398  tr * ( b * math::Pow2(sin(Py)) * c ),
399  tr * ( b * a * math::Pow2(sin(Pz)) ));
400 }
401 
402 
404 
405 
406 template<typename GridT, typename FieldT, typename InterruptT>
407 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
408  math::TemporalIntegrationScheme TemporalScheme>
409 inline
413  mParent(parent),
414  mVec(NULL),
415  mMinAbsV(1e-6),
416  mMap(*(parent.mTracker.grid().transform().template constMap<MapT>())),
417  mTask(0),
418  mIsMaster(true)
419 {
420 }
421 
422 template<typename GridT, typename FieldT, typename InterruptT>
423 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
424  math::TemporalIntegrationScheme TemporalScheme>
425 inline
426 LevelSetAdvection<GridT, FieldT, InterruptT>::
427 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
428 LevelSetAdvect(const LevelSetAdvect& other):
429  mParent(other.mParent),
430  mVec(other.mVec),
431  mMinAbsV(other.mMinAbsV),
432  mMaxAbsV(other.mMaxAbsV),
433  mMap(other.mMap),
434  mTask(other.mTask),
435  mIsMaster(false)
436 {
437 }
438 
439 template<typename GridT, typename FieldT, typename InterruptT>
440 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
441  math::TemporalIntegrationScheme TemporalScheme>
442 inline
443 LevelSetAdvection<GridT, FieldT, InterruptT>::
444 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
445 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
446  mParent(other.mParent),
447  mVec(other.mVec),
448  mMinAbsV(other.mMinAbsV),
449  mMaxAbsV(other.mMaxAbsV),
450  mMap(other.mMap),
451  mTask(other.mTask),
452  mIsMaster(false)
453 {
454 }
455 
456 template<typename GridT, typename FieldT, typename InterruptT>
457 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
458  math::TemporalIntegrationScheme TemporalScheme>
459 inline size_t
462 advect(ScalarType time0, ScalarType time1)
463 {
464  size_t count = 0;
465  while (time1 > time0 && mParent.mTracker.checkInterrupter()) {
467  mParent.mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme == math::TVD_RK3 ? 2 : 1);
468 
469  const ScalarType dt = this->sampleField(time0, time1);
470  if (math::isExactlyEqual(dt, ScalarType(0.0))) break;//V is essentially zero so terminate
471 
472  switch(TemporalScheme) {//switch is resolved at compile-time
473  case math::TVD_RK1:
474  // Perform one explicit Euler step: t1 = t0 + dt
475  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
476  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
477  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
478  this->cook(PARALLEL_FOR, 1);
479  break;
480  case math::TVD_RK2:
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 
487  // Convex combine explict Euler step: t2 = t0 + dt
488  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * V.Grad_t1(0))
489  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), /*phi=*/1, /*result=*/1);
490  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
491  this->cook(PARALLEL_FOR, 1);
492  break;
493  case math::TVD_RK3:
494  // Perform one explicit Euler step: t1 = t0 + dt
495  // Phi_t1(1) = Phi_t0(0) - dt * VdotG_t0(0)
496  mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, /*result=*/1);
497  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
498  this->cook(PARALLEL_FOR, 1);
499 
500  // Convex combine explict Euler step: t2 = t0 + dt/2
501  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * V.Grad_t1(0))
502  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), /*phi=*/1, /*result=*/2);
503  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
504  this->cook(PARALLEL_FOR, 2);
505 
506  // Convex combine explict Euler step: t3 = t0 + dt
507  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * V.Grad_t2(0)
508  mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), /*phi=*/1, /*result=*/2);
509  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
510  this->cook(PARALLEL_FOR, 2);
511  break;
512  default:
513  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
514  }
515 
516  time0 += dt;
517  ++count;
518  mParent.mTracker.mLeafs->removeAuxBuffers();
519  this->clearField();
520 
522  mParent.mTracker.track();
523  }
524  return count;
525 }
526 
527 template<typename GridT, typename FieldT, typename InterruptT>
528 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
529  math::TemporalIntegrationScheme TemporalScheme>
530 inline typename GridT::ValueType
531 LevelSetAdvection<GridT, FieldT, InterruptT>::
532 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
533 sampleField(ScalarType time0, ScalarType time1)
534 {
535  mMaxAbsV = mMinAbsV;
536  if (mParent.mTracker.mLeafs->leafCount()==0) return ScalarType(0.0);
537  mVec = new VectorType*[mParent.mTracker.mLeafs->leafCount()];
538  if (mParent.mField.transform() == mParent.mTracker.mGrid.transform()) {
539  mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0);
540  } else {
541  mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0);
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) : ScalarType(1.0))/math::Sqrt(ScalarType(3.0));
547  return math::Min(time1-time0, ScalarType(CFL*mParent.mTracker.voxelSize()/math::Sqrt(mMaxAbsV)));
548 }
549 
550 template<typename GridT, typename FieldT, typename InterruptT>
551 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
552  math::TemporalIntegrationScheme TemporalScheme>
553 inline void
554 LevelSetAdvection<GridT, FieldT, InterruptT>::
555 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
556 sampleXformedField(const RangeType& range, ScalarType time)
557 {
558  typedef typename LeafType::ValueOnCIter VoxelIterT;
559  const MapT& map = mMap;
560  mParent.mTracker.checkInterrupter();
561  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
562  const LeafType& leaf = mParent.mTracker.mLeafs->leaf(n);
563  VectorType* vec = new VectorType[leaf.onVoxelCount()];
564  int m = 0;
565  for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
566  const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time);
567  mMaxAbsV = math::Max(mMaxAbsV, ScalarType(math::Pow2(V[0])+math::Pow2(V[1])+math::Pow2(V[2])));
568  vec[m] = V;
569  }
570  mVec[n] = vec;
571  }
572 }
573 
574 template<typename GridT, typename FieldT, typename InterruptT>
575 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
576  math::TemporalIntegrationScheme TemporalScheme>
577 inline void
578 LevelSetAdvection<GridT, FieldT, InterruptT>::
579 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
580 sampleAlignedField(const RangeType& range, ScalarType time)
581 {
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.mLeafs->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(), time);
590  mMaxAbsV = math::Max(mMaxAbsV, ScalarType(math::Pow2(V[0])+math::Pow2(V[1])+math::Pow2(V[2])));
591  vec[m] = 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.mLeafs->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.mLeafs->getRange());
623  } else if (mode == PARALLEL_FOR) {
624  tbb::parallel_for(mParent.mTracker.mLeafs->getRange(mParent.mTracker.getGrainSize()), *this);
625  } else if (mode == PARALLEL_REDUCE) {
626  tbb::parallel_reduce(mParent.mTracker.mLeafs->getRange(mParent.mTracker.getGrainSize()), *this);
627  } else {
628  throw std::runtime_error("Undefined threading mode");
629  }
630 
631  mParent.mTracker.mLeafs->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  Stencil stencil(mParent.mTracker.mGrid);
652  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
653  BufferType& result = mParent.mTracker.mLeafs->getBuffer(n, resultBuffer);
654  const VectorType* vec = mVec[n];
655  int m=0;
656  for (VoxelIterT iter = mParent.mTracker.mLeafs->leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
657  stencil.moveTo(iter);
658  const VectorType V = vec[m], G = math::GradientBiased<MapT, SpatialScheme>::result(map, stencil, V);
659  result.setValue(iter.pos(), *iter - dt * V.dot(G));
660  }
661  }
662 }
663 
664 // Convex combination of Phi and a forward Euler advection steps:
665 // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * V.Grad(0));
666 template<typename GridT, typename FieldT, typename InterruptT>
667 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
668  math::TemporalIntegrationScheme TemporalScheme>
669 inline void
670 LevelSetAdvection<GridT, FieldT, InterruptT>::
671 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
672 euler2(const RangeType& range, ScalarType dt, ScalarType alpha, Index phiBuffer, Index resultBuffer)
673 {
674  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
675  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
676  typedef typename LeafType::ValueOnCIter VoxelIterT;
677  mParent.mTracker.checkInterrupter();
678  const MapT& map = mMap;
679  const ScalarType beta = ScalarType(1.0) - alpha;
680  Stencil stencil(mParent.mTracker.mGrid);
681  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
682  const BufferType& phi = mParent.mTracker.mLeafs->getBuffer(n, phiBuffer);
683  BufferType& result = mParent.mTracker.mLeafs->getBuffer(n, resultBuffer);
684  const VectorType* vec = mVec[n];
685  int m=0;
686  for (VoxelIterT iter = mParent.mTracker.mLeafs->leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
687  stencil.moveTo(iter);
688  const VectorType V = vec[m], G = math::GradientBiased<MapT, SpatialScheme>::result(map, stencil, V);
689  result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
690  }
691  }
692 }
693 
694 } // namespace tools
695 } // namespace OPENVDB_VERSION_NAME
696 } // namespace openvdb
697 
698 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
699 
700 // Copyright (c) 2012 DreamWorks Animation LLC
701 // All rights reserved. This software is distributed under the
702 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )