OpenVDB  7.0.0
LevelSetMorph.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
10 
11 #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
12 #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
13 
14 #include "LevelSetTracker.h"
15 #include "Interpolation.h" // for BoxSampler, etc.
17 #include <functional>
18 #include <limits>
19 
20 
21 namespace openvdb {
23 namespace OPENVDB_VERSION_NAME {
24 namespace tools {
25 
45 template<typename GridT,
46  typename InterruptT = util::NullInterrupter>
48 {
49 public:
50  using GridType = GridT;
51  using TreeType = typename GridT::TreeType;
53  using LeafRange = typename TrackerT::LeafRange;
54  using LeafType = typename TrackerT::LeafType;
55  using BufferType = typename TrackerT::BufferType;
56  using ValueType = typename TrackerT::ValueType;
57 
59  LevelSetMorphing(GridT& sourceGrid, const GridT& targetGrid, InterruptT* interrupt = nullptr)
60  : mTracker(sourceGrid, interrupt)
61  , mTarget(&targetGrid)
62  , mMask(nullptr)
63  , mSpatialScheme(math::HJWENO5_BIAS)
64  , mTemporalScheme(math::TVD_RK2)
65  , mMinMask(0)
66  , mDeltaMask(1)
67  , mInvertMask(false)
68  {
69  }
70 
71  virtual ~LevelSetMorphing() {}
72 
74  void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; }
75 
77  void setAlphaMask(const GridT& maskGrid) { mMask = &maskGrid; }
78 
80  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
82  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
83 
85  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
87  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
88 
91  {
92  return mTracker.getSpatialScheme();
93  }
96  {
97  mTracker.setSpatialScheme(scheme);
98  }
101  {
102  return mTracker.getTemporalScheme();
103  }
106  {
107  mTracker.setTemporalScheme(scheme);
108  }
110  int getNormCount() const { return mTracker.getNormCount(); }
112  void setNormCount(int n) { mTracker.setNormCount(n); }
113 
115  int getGrainSize() const { return mTracker.getGrainSize(); }
118  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
119 
122  ValueType minMask() const { return mMinMask; }
123 
126  ValueType maxMask() const { return mDeltaMask + mMinMask; }
127 
136  {
137  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
138  mMinMask = min;
139  mDeltaMask = max-min;
140  }
141 
144  bool isMaskInverted() const { return mInvertMask; }
147  void invertMask(bool invert=true) { mInvertMask = invert; }
148 
153  size_t advect(ValueType time0, ValueType time1);
154 
155 private:
156 
157  // disallow copy construction and copy by assignment!
158  LevelSetMorphing(const LevelSetMorphing&);// not implemented
159  LevelSetMorphing& operator=(const LevelSetMorphing&);// not implemented
160 
161  template<math::BiasedGradientScheme SpatialScheme>
162  size_t advect1(ValueType time0, ValueType time1);
163 
164  template<math::BiasedGradientScheme SpatialScheme,
165  math::TemporalIntegrationScheme TemporalScheme>
166  size_t advect2(ValueType time0, ValueType time1);
167 
168  template<math::BiasedGradientScheme SpatialScheme,
169  math::TemporalIntegrationScheme TemporalScheme,
170  typename MapType>
171  size_t advect3(ValueType time0, ValueType time1);
172 
173  TrackerT mTracker;
174  const GridT *mTarget, *mMask;
175  math::BiasedGradientScheme mSpatialScheme;
176  math::TemporalIntegrationScheme mTemporalScheme;
177  ValueType mMinMask, mDeltaMask;
178  bool mInvertMask;
179 
180  // This templated private class implements all the level set magic.
181  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
182  math::TemporalIntegrationScheme TemporalScheme>
183  struct Morph
184  {
188  Morph(const Morph& other);
190  Morph(Morph& other, tbb::split);
192  virtual ~Morph() {}
195  size_t advect(ValueType time0, ValueType time1);
197  void operator()(const LeafRange& r) const
198  {
199  if (mTask) mTask(const_cast<Morph*>(this), r);
200  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
201  }
203  void operator()(const LeafRange& r)
204  {
205  if (mTask) mTask(this, r);
206  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
207  }
209  void join(const Morph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
210 
212  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
213  // method calling tbb
214  void cook(ThreadingMode mode, size_t swapBuffer = 0);
215 
217  typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer);
218  void sampleXformedSpeed(const LeafRange& r, Index speedBuffer);
219  void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer);
220 
221  // Convex combination of Phi and a forward Euler advection steps:
222  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|);
223  template <int Nominator, int Denominator>
224  void euler(const LeafRange&, ValueType, Index, Index, Index);
225  inline void euler01(const LeafRange& r, ValueType t, Index s) {this->euler<0,1>(r,t,0,1,s);}
226  inline void euler12(const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
227  inline void euler34(const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
228  inline void euler13(const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
229 
230  using FuncType = typename std::function<void (Morph*, const LeafRange&)>;
231  LevelSetMorphing* mParent;
232  ValueType mMinAbsS, mMaxAbsS;
233  const MapT* mMap;
234  FuncType mTask;
235  }; // end of private Morph struct
236 
237 };//end of LevelSetMorphing
238 
239 template<typename GridT, typename InterruptT>
240 inline size_t
242 {
243  switch (mSpatialScheme) {
244  case math::FIRST_BIAS:
245  return this->advect1<math::FIRST_BIAS >(time0, time1);
246  //case math::SECOND_BIAS:
247  //return this->advect1<math::SECOND_BIAS >(time0, time1);
248  //case math::THIRD_BIAS:
249  //return this->advect1<math::THIRD_BIAS >(time0, time1);
250  //case math::WENO5_BIAS:
251  //return this->advect1<math::WENO5_BIAS >(time0, time1);
252  case math::HJWENO5_BIAS:
253  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
254  case math::SECOND_BIAS:
255  case math::THIRD_BIAS:
256  case math::WENO5_BIAS:
257  case math::UNKNOWN_BIAS:
258  default:
259  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
260  }
261  return 0;
262 }
263 
264 template<typename GridT, typename InterruptT>
265 template<math::BiasedGradientScheme SpatialScheme>
266 inline size_t
268 {
269  switch (mTemporalScheme) {
270  case math::TVD_RK1:
271  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
272  case math::TVD_RK2:
273  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
274  case math::TVD_RK3:
275  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
276  case math::UNKNOWN_TIS:
277  default:
278  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
279  }
280  return 0;
281 }
282 
283 template<typename GridT, typename InterruptT>
284 template<math::BiasedGradientScheme SpatialScheme,
285  math::TemporalIntegrationScheme TemporalScheme>
286 inline size_t
288 {
289  const math::Transform& trans = mTracker.grid().transform();
290  if (trans.mapType() == math::UniformScaleMap::mapType()) {
291  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
292  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
293  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
294  time0, time1);
295  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
296  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
297  } else if (trans.mapType() == math::TranslationMap::mapType()) {
298  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
299  } else {
300  OPENVDB_THROW(ValueError, "MapType not supported!");
301  }
302  return 0;
303 }
304 
305 template<typename GridT, typename InterruptT>
306 template<math::BiasedGradientScheme SpatialScheme,
307  math::TemporalIntegrationScheme TemporalScheme,
308  typename MapT>
309 inline size_t
311 {
312  Morph<MapT, SpatialScheme, TemporalScheme> tmp(*this);
313  return tmp.advect(time0, time1);
314 }
315 
316 
318 
319 template<typename GridT, typename InterruptT>
320 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
321  math::TemporalIntegrationScheme TemporalScheme>
322 inline
326  : mParent(&parent)
327  , mMinAbsS(ValueType(1e-6))
328  , mMap(parent.mTracker.grid().transform().template constMap<MapT>().get())
329  , mTask(nullptr)
330 {
331 }
332 
333 template<typename GridT, typename InterruptT>
334 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
335  math::TemporalIntegrationScheme TemporalScheme>
336 inline
339 Morph(const Morph& other)
340  : mParent(other.mParent)
341  , mMinAbsS(other.mMinAbsS)
342  , mMaxAbsS(other.mMaxAbsS)
343  , mMap(other.mMap)
344  , mTask(other.mTask)
345 {
346 }
347 
348 template<typename GridT, typename InterruptT>
349 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
350  math::TemporalIntegrationScheme TemporalScheme>
351 inline
354 Morph(Morph& other, tbb::split)
355  : mParent(other.mParent)
356  , mMinAbsS(other.mMinAbsS)
357  , mMaxAbsS(other.mMaxAbsS)
358  , mMap(other.mMap)
359  , mTask(other.mTask)
360 {
361 }
362 
363 template<typename GridT, typename InterruptT>
364 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
365  math::TemporalIntegrationScheme TemporalScheme>
366 inline size_t
369 advect(ValueType time0, ValueType time1)
370 {
371  namespace ph = std::placeholders;
372 
373  // Make sure we have enough temporal auxiliary buffers for the time
374  // integration AS WELL AS an extra buffer with the speed function!
375  static const Index auxBuffers = 1 + (TemporalScheme == math::TVD_RK3 ? 2 : 1);
376  size_t countCFL = 0;
377  while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
378  mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
379 
380  const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
381  if ( math::isZero(dt) ) break;//V is essentially zero so terminate
382 
383  OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time
384  switch(TemporalScheme) {
385  case math::TVD_RK1:
386  // Perform one explicit Euler step: t1 = t0 + dt
387  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
388  mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2);
389 
390  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
391  this->cook(PARALLEL_FOR, 1);
392  break;
393  case math::TVD_RK2:
394  // Perform one explicit Euler step: t1 = t0 + dt
395  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
396  mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2);
397 
398  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
399  this->cook(PARALLEL_FOR, 1);
400 
401  // Convex combine explict Euler step: t2 = t0 + dt
402  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|)
403  mTask = std::bind(&Morph::euler12, ph::_1, ph::_2, dt);
404 
405  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
406  this->cook(PARALLEL_FOR, 1);
407  break;
408  case math::TVD_RK3:
409  // Perform one explicit Euler step: t1 = t0 + dt
410  // Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]|
411  mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/3);
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 
416  // Convex combine explict Euler step: t2 = t0 + dt/2
417  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|)
418  mTask = std::bind(&Morph::euler34, ph::_1, ph::_2, dt);
419 
420  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
421  this->cook(PARALLEL_FOR, 2);
422 
423  // Convex combine explict Euler step: t3 = t0 + dt
424  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|)
425  mTask = std::bind(&Morph::euler13, ph::_1, ph::_2, dt);
426 
427  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
428  this->cook(PARALLEL_FOR, 2);
429  break;
430  case math::UNKNOWN_TIS:
431  default:
432  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
433  }//end of compile-time resolved switch
435 
436  time0 += dt;
437  ++countCFL;
438  mParent->mTracker.leafs().removeAuxBuffers();
439 
440  // Track the narrow band
441  mParent->mTracker.track();
442  }//end wile-loop over time
443 
444  return countCFL;//number of CLF propagation steps
445 }
446 
447 template<typename GridT, typename InterruptT>
448 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
449  math::TemporalIntegrationScheme TemporalScheme>
450 inline typename GridT::ValueType
453 sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer)
454 {
455  namespace ph = std::placeholders;
456 
457  mMaxAbsS = mMinAbsS;
458  const size_t leafCount = mParent->mTracker.leafs().leafCount();
459  if (leafCount==0 || time0 >= time1) return ValueType(0);
460 
461  const math::Transform& xform = mParent->mTracker.grid().transform();
462  if (mParent->mTarget->transform() == xform &&
463  (mParent->mMask == nullptr || mParent->mMask->transform() == xform)) {
464  mTask = std::bind(&Morph::sampleAlignedSpeed, ph::_1, ph::_2, speedBuffer);
465  } else {
466  mTask = std::bind(&Morph::sampleXformedSpeed, ph::_1, ph::_2, speedBuffer);
467  }
468  this->cook(PARALLEL_REDUCE);
469  if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ValueType(0);//speed is essentially zero
470  static const ValueType CFL = (TemporalScheme == math::TVD_RK1 ? ValueType(0.3) :
471  TemporalScheme == math::TVD_RK2 ? ValueType(0.9) :
472  ValueType(1.0))/math::Sqrt(ValueType(3.0));
473  const ValueType dt = math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
474  return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
475 }
476 
477 template<typename GridT, typename InterruptT>
478 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
479  math::TemporalIntegrationScheme TemporalScheme>
480 inline void
483 sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
484 {
485  using VoxelIterT = typename LeafType::ValueOnCIter;
487 
488  const MapT& map = *mMap;
489  mParent->mTracker.checkInterrupter();
490 
491  typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
492  SamplerT target(targetAcc, mParent->mTarget->transform());
493  if (mParent->mMask == nullptr) {
494  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
495  ValueType* speed = leafIter.buffer(speedBuffer).data();
496  bool isZero = true;
497  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
498  ValueType& s = speed[voxelIter.pos()];
499  s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
500  if (!math::isApproxZero(s)) isZero = false;
501  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
502  }
503  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
504  }
505  } else {
506  const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
507  const bool invMask = mParent->isMaskInverted();
508  typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
509  SamplerT mask(maskAcc, mParent->mMask->transform());
510  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
511  ValueType* speed = leafIter.buffer(speedBuffer).data();
512  bool isZero = true;
513  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
514  const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());//world space
515  const ValueType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm);
516  ValueType& s = speed[voxelIter.pos()];
517  s -= target.wsSample(xyz);
518  s *= invMask ? 1 - a : a;
519  if (!math::isApproxZero(s)) isZero = false;
520  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
521  }
522  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
523  }
524  }
525 }
526 
527 template<typename GridT, typename InterruptT>
528 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
529  math::TemporalIntegrationScheme TemporalScheme>
530 inline void
533 sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
534 {
535  using VoxelIterT = typename LeafType::ValueOnCIter;
536 
537  mParent->mTracker.checkInterrupter();
538 
539  typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
540 
541  if (mParent->mMask == nullptr) {
542  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
543  ValueType* speed = leafIter.buffer(speedBuffer).data();
544  bool isZero = true;
545  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
546  ValueType& s = speed[voxelIter.pos()];
547  s -= target.getValue(voxelIter.getCoord());
548  if (!math::isApproxZero(s)) isZero = false;
549  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
550  }
551  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
552  }
553  } else {
554  const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
555  const bool invMask = mParent->isMaskInverted();
556  typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
557  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
558  ValueType* speed = leafIter.buffer(speedBuffer).data();
559  bool isZero = true;
560  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
561  const Coord ijk = voxelIter.getCoord();//index space
562  const ValueType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm);
563  ValueType& s = speed[voxelIter.pos()];
564  s -= target.getValue(ijk);
565  s *= invMask ? 1 - a : a;
566  if (!math::isApproxZero(s)) isZero = false;
567  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
568  }
569  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
570  }
571  }
572 }
573 
574 template<typename GridT, typename InterruptT>
575 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
576  math::TemporalIntegrationScheme TemporalScheme>
577 inline void
580 cook(ThreadingMode mode, size_t swapBuffer)
581 {
582  mParent->mTracker.startInterrupter("Morphing level set");
583 
584  const int grainSize = mParent->mTracker.getGrainSize();
585  const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
586 
587  if (mParent->mTracker.getGrainSize()==0) {
588  (*this)(range);
589  } else if (mode == PARALLEL_FOR) {
590  tbb::parallel_for(range, *this);
591  } else if (mode == PARALLEL_REDUCE) {
592  tbb::parallel_reduce(range, *this);
593  } else {
594  OPENVDB_THROW(ValueError, "expected threading mode " << int(PARALLEL_FOR)
595  << " or " << int(PARALLEL_REDUCE) << ", got " << int(mode));
596  }
597 
598  mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
599 
600  mParent->mTracker.endInterrupter();
601 }
602 
603 template<typename GridT, typename InterruptT>
604 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
605  math::TemporalIntegrationScheme TemporalScheme>
606 template <int Nominator, int Denominator>
607 inline void
610 euler(const LeafRange& range, ValueType dt,
611  Index phiBuffer, Index resultBuffer, Index speedBuffer)
612 {
613  using SchemeT = math::BIAS_SCHEME<SpatialScheme>;
614  using StencilT = typename SchemeT::template ISStencil<GridType>::StencilType;
615  using VoxelIterT = typename LeafType::ValueOnCIter;
617 
618  static const ValueType Alpha = ValueType(Nominator)/ValueType(Denominator);
619  static const ValueType Beta = ValueType(1) - Alpha;
620 
621  mParent->mTracker.checkInterrupter();
622  const MapT& map = *mMap;
623  StencilT stencil(mParent->mTracker.grid());
624 
625  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
626  const ValueType* speed = leafIter.buffer(speedBuffer).data();
628  const ValueType* phi = leafIter.buffer(phiBuffer).data();
629  ValueType* result = leafIter.buffer(resultBuffer).data();
630  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
631  const Index n = voxelIter.pos();
632  if (math::isApproxZero(speed[n])) continue;
633  stencil.moveTo(voxelIter);
634  const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
635  result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
636  }//loop over active voxels in the leaf of the mask
637  }//loop over leafs of the level set
638 }
639 
640 } // namespace tools
641 } // namespace OPENVDB_VERSION_NAME
642 } // namespace openvdb
643 
644 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:281
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
Definition: FiniteDifference.h:236
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:351
typename TrackerT::LeafType LeafType
Definition: LevelSetMorph.h:54
typename TrackerT::LeafRange LeafRange
Definition: LevelSetMorph.h:53
ValueType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetMorph.h:122
bool checkInterrupter()
Definition: LevelSetTracker.h:382
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:144
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
Definition: FiniteDifference.h:238
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:234
void setTarget(const GridT &targetGrid)
Redefine the target level set.
Definition: LevelSetMorph.h:74
GridT GridType
Definition: LevelSetMorph.h:50
int getNormCount() const
Return the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:110
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=nullptr)
Main constructor.
Definition: LevelSetMorph.h:59
typename GridT::TreeType TreeType
Definition: LevelSetMorph.h:51
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition: LevelSetMorph.h:47
Definition: FiniteDifference.h:237
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:542
typename TrackerT::BufferType BufferType
Definition: LevelSetMorph.h:55
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:603
Definition: FiniteDifference.h:168
Definition: Exceptions.h:65
Definition: FiniteDifference.h:235
ValueType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetMorph.h:126
Definition: FiniteDifference.h:166
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:112
Definition: FiniteDifference.h:169
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:95
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
Definition: Operators.h:854
virtual ~LevelSetMorphing()
Definition: LevelSetMorph.h:71
typename TreeType::ValueType ValueType
Definition: LevelSetTracker.h:63
Definition: Operators.h:126
int getGrainSize() const
Return the grain size used for multithreading.
Definition: LevelSetMorph.h:115
void startInterrupter(const char *msg)
Definition: LevelSetTracker.h:366
Definition: Exceptions.h:13
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
SIMD Intrinsic Headers.
Definition: Platform.h:114
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:55
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:708
Definition: FiniteDifference.h:170
Definition: FiniteDifference.h:171
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:82
Definition: Transform.h:39
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:85
Index32 Index
Definition: Types.h:31
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:115
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x².
Definition: Math.h:229
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:87
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:100
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:80
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:102
typename TreeType::LeafNodeType LeafType
Definition: LevelSetTracker.h:62
typename TrackerT::ValueType ValueType
Definition: LevelSetMorph.h:56
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:293
Definition: FiniteDifference.h:167
typename LeafManagerType::LeafRange LeafRange
Definition: LevelSetTracker.h:65
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:388
Name mapType() const
Return the transformation map&#39;s type-name.
Definition: Transform.h:66
Class that provides the interface for continuous sampling of values in a tree.
Definition: Interpolation.h:283
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:90
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
void setAlphaMask(const GridT &maskGrid)
Define the alpha mask.
Definition: LevelSetMorph.h:77
typename LeafManagerType::BufferType BufferType
Definition: LevelSetTracker.h:66
void setGrainSize(int grainsize)
Set the grain size used for multithreading.
Definition: LevelSetMorph.h:118
void setMaskRange(ValueType min, ValueType max)
Define the range for the (optional) scalar mask.
Definition: LevelSetMorph.h:135
LeafManagerType & leafs()
Definition: LevelSetTracker.h:187
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:105
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:165
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:147