OpenVDB  4.0.1
LevelSetMorph.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2017 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 
37 
38 #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
39 #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
40 
41 #include "LevelSetTracker.h"
42 #include "Interpolation.h" // for BoxSampler, etc.
44 
45 namespace openvdb {
47 namespace OPENVDB_VERSION_NAME {
48 namespace tools {
49 
50 
70 template<typename GridT,
71  typename InterruptT = util::NullInterrupter>
73 {
74 public:
75  typedef GridT GridType;
76  typedef typename GridT::TreeType TreeType;
78  typedef typename TrackerT::LeafRange LeafRange;
79  typedef typename TrackerT::LeafType LeafType;
80  typedef typename TrackerT::BufferType BufferType;
81  typedef typename TrackerT::ValueType ValueType;
82 
84  LevelSetMorphing(GridT& sourceGrid, const GridT& targetGrid, InterruptT* interrupt = nullptr)
85  : mTracker(sourceGrid, interrupt)
86  , mTarget(&targetGrid)
87  , mMask(nullptr)
88  , mSpatialScheme(math::HJWENO5_BIAS)
89  , mTemporalScheme(math::TVD_RK2)
90  , mMinMask(0)
91  , mDeltaMask(1)
92  , mInvertMask(false)
93  {
94  }
95 
96  virtual ~LevelSetMorphing() {}
97 
99  void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; }
100 
102  void setAlphaMask(const GridT& maskGrid) { mMask = &maskGrid; }
103 
105  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
107  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
108 
110  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
112  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
113 
116  {
117  return mTracker.getSpatialScheme();
118  }
121  {
122  mTracker.setSpatialScheme(scheme);
123  }
126  {
127  return mTracker.getTemporalScheme();
128  }
131  {
132  mTracker.setTemporalScheme(scheme);
133  }
135  int getNormCount() const { return mTracker.getNormCount(); }
137  void setNormCount(int n) { mTracker.setNormCount(n); }
138 
140  int getGrainSize() const { return mTracker.getGrainSize(); }
143  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
144 
147  ValueType minMask() const { return mMinMask; }
148 
151  ValueType maxMask() const { return mDeltaMask + mMinMask; }
152 
160  void setMaskRange(ValueType min, ValueType max)
161  {
162  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
163  mMinMask = min;
164  mDeltaMask = max-min;
165  }
166 
169  bool isMaskInverted() const { return mInvertMask; }
172  void invertMask(bool invert=true) { mInvertMask = invert; }
173 
178  size_t advect(ValueType time0, ValueType time1);
179 
180 private:
181 
182  // disallow copy construction and copy by assignment!
183  LevelSetMorphing(const LevelSetMorphing&);// not implemented
184  LevelSetMorphing& operator=(const LevelSetMorphing&);// not implemented
185 
186  template<math::BiasedGradientScheme SpatialScheme>
187  size_t advect1(ValueType time0, ValueType time1);
188 
189  template<math::BiasedGradientScheme SpatialScheme,
190  math::TemporalIntegrationScheme TemporalScheme>
191  size_t advect2(ValueType time0, ValueType time1);
192 
193  template<math::BiasedGradientScheme SpatialScheme,
194  math::TemporalIntegrationScheme TemporalScheme,
195  typename MapType>
196  size_t advect3(ValueType time0, ValueType time1);
197 
198  TrackerT mTracker;
199  const GridT *mTarget, *mMask;
200  math::BiasedGradientScheme mSpatialScheme;
201  math::TemporalIntegrationScheme mTemporalScheme;
202  ValueType mMinMask, mDeltaMask;
203  bool mInvertMask;
204 
205  // This templated private class implements all the level set magic.
206  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
207  math::TemporalIntegrationScheme TemporalScheme>
208  struct Morph
209  {
213  Morph(const Morph& other);
215  Morph(Morph& other, tbb::split);
217  virtual ~Morph() {}
220  size_t advect(ValueType time0, ValueType time1);
222  void operator()(const LeafRange& r) const
223  {
224  if (mTask) mTask(const_cast<Morph*>(this), r);
225  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
226  }
228  void operator()(const LeafRange& r)
229  {
230  if (mTask) mTask(this, r);
231  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
232  }
234  void join(const Morph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
235 
237  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
238  // method calling tbb
239  void cook(ThreadingMode mode, size_t swapBuffer = 0);
240 
242  typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer);
243  void sampleXformedSpeed(const LeafRange& r, Index speedBuffer);
244  void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer);
245 
246  // Convex combination of Phi and a forward Euler advection steps:
247  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|);
248  template <int Nominator, int Denominator>
249  void euler(const LeafRange&, ValueType, Index, Index, Index);
250  inline void euler01(const LeafRange& r, ValueType t, Index s) {this->euler<0,1>(r,t,0,1,s);}
251  inline void euler12(const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
252  inline void euler34(const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
253  inline void euler13(const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
254 
255  typedef typename boost::function<void (Morph*, const LeafRange&)> FuncType;
256  LevelSetMorphing* mParent;
257  ValueType mMinAbsS, mMaxAbsS;
258  const MapT* mMap;
259  FuncType mTask;
260  }; // end of private Morph struct
261 
262 };//end of LevelSetMorphing
263 
264 template<typename GridT, typename InterruptT>
265 inline size_t
267 {
268  switch (mSpatialScheme) {
269  case math::FIRST_BIAS:
270  return this->advect1<math::FIRST_BIAS >(time0, time1);
271  //case math::SECOND_BIAS:
272  //return this->advect1<math::SECOND_BIAS >(time0, time1);
273  //case math::THIRD_BIAS:
274  //return this->advect1<math::THIRD_BIAS >(time0, time1);
275  //case math::WENO5_BIAS:
276  //return this->advect1<math::WENO5_BIAS >(time0, time1);
277  case math::HJWENO5_BIAS:
278  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
279  case math::SECOND_BIAS:
280  case math::THIRD_BIAS:
281  case math::WENO5_BIAS:
282  case math::UNKNOWN_BIAS:
283  default:
284  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
285  }
286  return 0;
287 }
288 
289 template<typename GridT, typename InterruptT>
290 template<math::BiasedGradientScheme SpatialScheme>
291 inline size_t
293 {
294  switch (mTemporalScheme) {
295  case math::TVD_RK1:
296  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
297  case math::TVD_RK2:
298  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
299  case math::TVD_RK3:
300  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
301  case math::UNKNOWN_TIS:
302  default:
303  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
304  }
305  return 0;
306 }
307 
308 template<typename GridT, typename InterruptT>
309 template<math::BiasedGradientScheme SpatialScheme,
310  math::TemporalIntegrationScheme TemporalScheme>
311 inline size_t
313 {
314  const math::Transform& trans = mTracker.grid().transform();
315  if (trans.mapType() == math::UniformScaleMap::mapType()) {
316  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
317  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
318  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
319  time0, time1);
320  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
321  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
322  } else if (trans.mapType() == math::TranslationMap::mapType()) {
323  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
324  } else {
325  OPENVDB_THROW(ValueError, "MapType not supported!");
326  }
327  return 0;
328 }
329 
330 template<typename GridT, typename InterruptT>
331 template<math::BiasedGradientScheme SpatialScheme,
332  math::TemporalIntegrationScheme TemporalScheme,
333  typename MapT>
334 inline size_t
336 {
337  Morph<MapT, SpatialScheme, TemporalScheme> tmp(*this);
338  return tmp.advect(time0, time1);
339 }
340 
341 
343 
344 template<typename GridT, typename InterruptT>
345 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
346  math::TemporalIntegrationScheme TemporalScheme>
347 inline
351  : mParent(&parent)
352  , mMinAbsS(ValueType(1e-6))
353  , mMap(parent.mTracker.grid().transform().template constMap<MapT>().get())
354  , mTask(0)
355 {
356 }
357 
358 template<typename GridT, typename InterruptT>
359 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
360  math::TemporalIntegrationScheme TemporalScheme>
361 inline
364 Morph(const Morph& other)
365  : mParent(other.mParent)
366  , mMinAbsS(other.mMinAbsS)
367  , mMaxAbsS(other.mMaxAbsS)
368  , mMap(other.mMap)
369  , mTask(other.mTask)
370 {
371 }
372 
373 template<typename GridT, typename InterruptT>
374 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
375  math::TemporalIntegrationScheme TemporalScheme>
376 inline
379 Morph(Morph& other, tbb::split)
380  : mParent(other.mParent)
381  , mMinAbsS(other.mMinAbsS)
382  , mMaxAbsS(other.mMaxAbsS)
383  , mMap(other.mMap)
384  , mTask(other.mTask)
385 {
386 }
387 
388 template<typename GridT, typename InterruptT>
389 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
390  math::TemporalIntegrationScheme TemporalScheme>
391 inline size_t
394 advect(ValueType time0, ValueType time1)
395 {
396  // Make sure we have enough temporal auxiliary buffers for the time
397  // integration AS WELL AS an extra buffer with the speed function!
398  static const Index auxBuffers = 1 + (TemporalScheme == math::TVD_RK3 ? 2 : 1);
399  size_t countCFL = 0;
400  while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
401  mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
402 
403  const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
404  if ( math::isZero(dt) ) break;//V is essentially zero so terminate
405 
406  OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time
407  switch(TemporalScheme) {
408  case math::TVD_RK1:
409  // Perform one explicit Euler step: t1 = t0 + dt
410  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
411  mTask = boost::bind(&Morph::euler01, _1, _2, dt, /*speed*/2);
412 
413  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
414  this->cook(PARALLEL_FOR, 1);
415  break;
416  case math::TVD_RK2:
417  // Perform one explicit Euler step: t1 = t0 + dt
418  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
419  mTask = boost::bind(&Morph::euler01, _1, _2, dt, /*speed*/2);
420 
421  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
422  this->cook(PARALLEL_FOR, 1);
423 
424  // Convex combine explict Euler step: t2 = t0 + dt
425  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|)
426  mTask = boost::bind(&Morph::euler12, _1, _2, dt);
427 
428  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
429  this->cook(PARALLEL_FOR, 1);
430  break;
431  case math::TVD_RK3:
432  // Perform one explicit Euler step: t1 = t0 + dt
433  // Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]|
434  mTask = boost::bind(&Morph::euler01, _1, _2, dt, /*speed*/3);
435 
436  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
437  this->cook(PARALLEL_FOR, 1);
438 
439  // Convex combine explict Euler step: t2 = t0 + dt/2
440  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|)
441  mTask = boost::bind(&Morph::euler34, _1, _2, dt);
442 
443  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
444  this->cook(PARALLEL_FOR, 2);
445 
446  // Convex combine explict Euler step: t3 = t0 + dt
447  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|)
448  mTask = boost::bind(&Morph::euler13, _1, _2, dt);
449 
450  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
451  this->cook(PARALLEL_FOR, 2);
452  break;
453  case math::UNKNOWN_TIS:
454  default:
455  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
456  }//end of compile-time resolved switch
458 
459  time0 += dt;
460  ++countCFL;
461  mParent->mTracker.leafs().removeAuxBuffers();
462 
463  // Track the narrow band
464  mParent->mTracker.track();
465  }//end wile-loop over time
466 
467  return countCFL;//number of CLF propagation steps
468 }
469 
470 template<typename GridT, typename InterruptT>
471 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
472  math::TemporalIntegrationScheme TemporalScheme>
473 inline typename GridT::ValueType
476 sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer)
477 {
478  mMaxAbsS = mMinAbsS;
479  const size_t leafCount = mParent->mTracker.leafs().leafCount();
480  if (leafCount==0 || time0 >= time1) return ValueType(0);
481 
482  const math::Transform& xform = mParent->mTracker.grid().transform();
483  if (mParent->mTarget->transform() == xform &&
484  (mParent->mMask == nullptr || mParent->mMask->transform() == xform)) {
485  mTask = boost::bind(&Morph::sampleAlignedSpeed, _1, _2, speedBuffer);
486  } else {
487  mTask = boost::bind(&Morph::sampleXformedSpeed, _1, _2, speedBuffer);
488  }
489  this->cook(PARALLEL_REDUCE);
490  if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ValueType(0);//speed is essentially zero
491  static const ValueType CFL = (TemporalScheme == math::TVD_RK1 ? ValueType(0.3) :
492  TemporalScheme == math::TVD_RK2 ? ValueType(0.9) :
493  ValueType(1.0))/math::Sqrt(ValueType(3.0));
494  const ValueType dt = math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
495  return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
496 }
497 
498 template<typename GridT, typename InterruptT>
499 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
500  math::TemporalIntegrationScheme TemporalScheme>
501 inline void
504 sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
505 {
506  typedef typename LeafType::ValueOnCIter VoxelIterT;
508  const MapT& map = *mMap;
509  mParent->mTracker.checkInterrupter();
510 
511  typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
512  SamplerT target(targetAcc, mParent->mTarget->transform());
513  if (mParent->mMask == nullptr) {
514  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
515  ValueType* speed = leafIter.buffer(speedBuffer).data();
516  bool isZero = true;
517  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
518  ValueType& s = speed[voxelIter.pos()];
519  s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
520  if (!math::isApproxZero(s)) isZero = false;
521  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
522  }
523  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
524  }
525  } else {
526  const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
527  const bool invMask = mParent->isMaskInverted();
528  typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
529  SamplerT mask(maskAcc, mParent->mMask->transform());
530  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
531  ValueType* speed = leafIter.buffer(speedBuffer).data();
532  bool isZero = true;
533  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
534  const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());//world space
535  const ValueType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm);
536  ValueType& s = speed[voxelIter.pos()];
537  s -= target.wsSample(xyz);
538  s *= invMask ? 1 - a : a;
539  if (!math::isApproxZero(s)) isZero = false;
540  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
541  }
542  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
543  }
544  }
545 }
546 
547 template<typename GridT, typename InterruptT>
548 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
549  math::TemporalIntegrationScheme TemporalScheme>
550 inline void
553 sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
554 {
555  typedef typename LeafType::ValueOnCIter VoxelIterT;
556  mParent->mTracker.checkInterrupter();
557 
558  typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
559 
560  if (mParent->mMask == nullptr) {
561  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
562  ValueType* speed = leafIter.buffer(speedBuffer).data();
563  bool isZero = true;
564  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
565  ValueType& s = speed[voxelIter.pos()];
566  s -= target.getValue(voxelIter.getCoord());
567  if (!math::isApproxZero(s)) isZero = false;
568  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
569  }
570  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
571  }
572  } else {
573  const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
574  const bool invMask = mParent->isMaskInverted();
575  typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
576  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
577  ValueType* speed = leafIter.buffer(speedBuffer).data();
578  bool isZero = true;
579  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
580  const Coord ijk = voxelIter.getCoord();//index space
581  const ValueType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm);
582  ValueType& s = speed[voxelIter.pos()];
583  s -= target.getValue(ijk);
584  s *= invMask ? 1 - a : a;
585  if (!math::isApproxZero(s)) isZero = false;
586  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
587  }
588  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
589  }
590  }
591 }
592 
593 template<typename GridT, typename InterruptT>
594 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
595  math::TemporalIntegrationScheme TemporalScheme>
596 inline void
599 cook(ThreadingMode mode, size_t swapBuffer)
600 {
601  mParent->mTracker.startInterrupter("Morphing level set");
602 
603  const int grainSize = mParent->mTracker.getGrainSize();
604  const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
605 
606  if (mParent->mTracker.getGrainSize()==0) {
607  (*this)(range);
608  } else if (mode == PARALLEL_FOR) {
609  tbb::parallel_for(range, *this);
610  } else if (mode == PARALLEL_REDUCE) {
611  tbb::parallel_reduce(range, *this);
612  } else {
613  throw std::runtime_error("Undefined threading mode");
614  }
615 
616  mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
617 
618  mParent->mTracker.endInterrupter();
619 }
620 
621 template<typename GridT, typename InterruptT>
622 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
623  math::TemporalIntegrationScheme TemporalScheme>
624 template <int Nominator, int Denominator>
625 inline void
628 euler(const LeafRange& range, ValueType dt,
629  Index phiBuffer, Index resultBuffer, Index speedBuffer)
630 {
631  typedef math::BIAS_SCHEME<SpatialScheme> SchemeT;
632  typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
633  typedef typename LeafType::ValueOnCIter VoxelIterT;
635 
636  static const ValueType Alpha = ValueType(Nominator)/ValueType(Denominator);
637  static const ValueType Beta = ValueType(1) - Alpha;
638 
639  mParent->mTracker.checkInterrupter();
640  const MapT& map = *mMap;
641  StencilT stencil(mParent->mTracker.grid());
642 
643  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
644  const ValueType* speed = leafIter.buffer(speedBuffer).data();
646  const ValueType* phi = leafIter.buffer(phiBuffer).data();
647  ValueType* result = leafIter.buffer(resultBuffer).data();
648  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
649  const Index n = voxelIter.pos();
650  if (math::isApproxZero(speed[n])) continue;
651  stencil.moveTo(voxelIter);
652  const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
653  result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
654  }//loop over active voxels in the leaf of the mask
655  }//loop over leafs of the level set
656 }
657 
658 } // namespace tools
659 } // namespace OPENVDB_VERSION_NAME
660 } // namespace openvdb
661 
662 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
663 
664 // Copyright (c) 2012-2017 DreamWorks Animation LLC
665 // All rights reserved. This software is distributed under the
666 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:107
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:115
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:370
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:333
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
TrackerT::BufferType BufferType
Definition: LevelSetMorph.h:80
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
Class that provides the interface for continuous sampling of values in a tree.
Definition: Interpolation.h:310
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:112
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:272
TrackerT::LeafType LeafType
Definition: LevelSetMorph.h:79
void setGrainSize(int grainsize)
Set the grain size used for multithreading.
Definition: LevelSetMorph.h:143
Coord Abs(const Coord &xyz)
Definition: Coord.h:254
Definition: FiniteDifference.h:262
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:622
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:101
bool isMaskInverted() const
Return true if the mask is inverted, i.e. min->max in the original mask maps to 1->0 in the inverted ...
Definition: LevelSetMorph.h:169
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition: LevelSetMorph.h:72
Definition: Operators.h:152
Name mapType() const
Return the transformation map&#39;s type-name.
Definition: Transform.h:93
TreeType::LeafNodeType LeafType
Definition: LevelSetTracker.h:72
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:48
Definition: FiniteDifference.h:195
GridT::TreeType TreeType
Definition: LevelSetMorph.h:76
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
LevelSetTracker< GridT, InterruptT > TrackerT
Definition: LevelSetMorph.h:77
void setAlphaMask(const GridT &maskGrid)
Define the alpha mask.
Definition: LevelSetMorph.h:102
void startInterrupter(const char *msg)
Definition: LevelSetTracker.h:359
void setTarget(const GridT &targetGrid)
Redefine the target level set.
Definition: LevelSetMorph.h:99
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:105
int getNormCount() const
Return the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:135
void setMaskRange(ValueType min, ValueType max)
Define the range for the (optional) scalar mask.
Definition: LevelSetMorph.h:160
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:727
TrackerT::LeafRange LeafRange
Definition: LevelSetMorph.h:78
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:128
#define OPENVDB_VERSION_NAME
Definition: version.h:43
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:137
void invertMask(bool invert=true)
Invert the optional mask, i.e. min->max in the original mask maps to 1->0 in the inverted alpha mask...
Definition: LevelSetMorph.h:172
Definition: FiniteDifference.h:264
Definition: Exceptions.h:39
ValueType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetMorph.h:147
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:129
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:125
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: FiniteDifference.h:194
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
Definition: FiniteDifference.h:263
Definition: FiniteDifference.h:196
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=nullptr)
Main constructor.
Definition: LevelSetMorph.h:84
void rebuildAuxBuffers(size_t auxBuffersPerLeaf, bool serial=false)
Change the number of auxiliary buffers.
Definition: LeafManager.h:312
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:110
Definition: Exceptions.h:92
Definition: FiniteDifference.h:265
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:130
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
Definition: FiniteDifference.h:198
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
Iterator begin() const
Definition: LeafManager.h:186
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:336
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:67
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:120
Definition: FiniteDifference.h:193
LeafManagerType & leafs()
Definition: LevelSetTracker.h:186
ValueType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetMorph.h:151
Definition: Operators.h:860
int getGrainSize() const
Return the grain size used for multithreading.
Definition: LevelSetMorph.h:140
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Index32 Index
Definition: Types.h:57
TrackerT::ValueType ValueType
Definition: LevelSetMorph.h:81
LeafManagerType::BufferType BufferType
Definition: LevelSetTracker.h:76
GridT GridType
Definition: LevelSetMorph.h:75
Definition: FiniteDifference.h:197
TreeType::ValueType ValueType
Definition: LevelSetTracker.h:73
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
bool checkInterrupter()
Definition: LevelSetTracker.h:375
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:130
virtual ~LevelSetMorphing()
Definition: LevelSetMorph.h:96