HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LevelSetMorph.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
4 /// @author Ken Museth
5 ///
6 /// @file tools/LevelSetMorph.h
7 ///
8 /// @brief Shape morphology of level sets. Morphing from a source
9 /// narrow-band level sets to a target narrow-band level set.
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 
26 /// @brief Shape morphology of level sets. Morphing from a source
27 /// narrow-band level sets to a target narrow-band level set.
28 ///
29 /// @details
30 /// The @c InterruptType template argument below refers to any class
31 /// with the following interface:
32 /// @code
33 /// class Interrupter {
34 /// ...
35 /// public:
36 /// void start(const char* name = nullptr) // called when computations begin
37 /// void end() // called when computations end
38 /// bool wasInterrupted(int percent=-1) // return true to break computation
39 /// };
40 /// @endcode
41 ///
42 /// @note If no template argument is provided for this InterruptType,
43 /// the util::NullInterrupter is used, which implies that all interrupter
44 /// calls are no-ops (i.e., they incur no computational overhead).
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 
58  /// Main constructor
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 
73  /// Redefine the target level set
74  void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; }
75 
76  /// Define the alpha mask
77  void setAlphaMask(const GridT& maskGrid) { mMask = &maskGrid; }
78 
79  /// Return the spatial finite-difference scheme
80  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
81  /// Set the spatial finite-difference scheme
82  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
83 
84  /// Return the temporal integration scheme
85  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
86  /// Set the temporal integration scheme
87  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
88 
89  /// Return the spatial finite-difference scheme
91  {
92  return mTracker.getSpatialScheme();
93  }
94  /// Set the spatial finite-difference scheme
96  {
97  mTracker.setSpatialScheme(scheme);
98  }
99  /// Return the temporal integration scheme
101  {
102  return mTracker.getTemporalScheme();
103  }
104  /// Set the temporal integration scheme
106  {
107  mTracker.setTemporalScheme(scheme);
108  }
109  /// Return the number of normalizations performed per track or normalize call.
110  int getNormCount() const { return mTracker.getNormCount(); }
111  /// Set the number of normalizations performed per track or normalize call.
112  void setNormCount(int n) { mTracker.setNormCount(n); }
113 
114  /// Return the grain size used for multithreading
115  int getGrainSize() const { return mTracker.getGrainSize(); }
116  /// @brief Set the grain size used for multithreading.
117  /// @note A grain size of 0 or less disables multithreading!
118  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
119 
120  /// @brief Return the minimum value of the mask to be used for the
121  /// derivation of a smooth alpha value.
122  ValueType minMask() const { return mMinMask; }
123 
124  /// @brief Return the maximum value of the mask to be used for the
125  /// derivation of a smooth alpha value.
126  ValueType maxMask() const { return mDeltaMask + mMinMask; }
127 
128  /// @brief Define the range for the (optional) scalar mask.
129  /// @param min Minimum value of the range.
130  /// @param max Maximum value of the range.
131  /// @details Mask values outside the range maps to alpha values of
132  /// respectfully zero and one, and values inside the range maps
133  /// smoothly to 0->1 (unless of course the mask is inverted).
134  /// @throw ValueError if @a min is not smaller than @a max.
136  {
137  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
138  mMinMask = min;
139  mDeltaMask = max-min;
140  }
141 
142  /// @brief Return true if the mask is inverted, i.e. min->max in the
143  /// original mask maps to 1->0 in the inverted alpha mask.
144  bool isMaskInverted() const { return mInvertMask; }
145  /// @brief Invert the optional mask, i.e. min->max in the original
146  /// mask maps to 1->0 in the inverted alpha mask.
147  void invertMask(bool invert=true) { mInvertMask = invert; }
148 
149  /// @brief Advect the level set from its current time, @a time0, to its
150  /// final time, @a time1. If @a time0 > @a time1, perform backward advection.
151  ///
152  /// @return the number of CFL iterations used to advect from @a time0 to @a time1
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  {
185  /// Main constructor
187  /// Shallow copy constructor called by tbb::parallel_for() threads
188  Morph(const Morph& other);
189  /// Shallow copy constructor called by tbb::parallel_reduce() threads
190  Morph(Morph& other, tbb::split);
191  /// destructor
192  virtual ~Morph() {}
193  /// Advect the level set from its current time, time0, to its final time, time1.
194  /// @return number of CFL iterations
195  size_t advect(ValueType time0, ValueType time1);
196  /// Used internally by tbb::parallel_for()
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  }
202  /// Used internally by tbb::parallel_reduce()
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  }
208  /// This is only called by tbb::parallel_reduce() threads
209  void join(const Morph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
210 
211  /// Enum to define the type of multithreading
212  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
213  // method calling tbb
214  void cook(ThreadingMode mode, size_t swapBuffer = 0);
215 
216  /// Sample field and return the CFT time step
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
287 LevelSetMorphing<GridT, InterruptT>::advect2(ValueType time0, ValueType time1)
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
310 LevelSetMorphing<GridT, InterruptT>::advect3(ValueType time0, ValueType time1)
311 {
312  Morph<MapT, SpatialScheme, TemporalScheme> tmp(*this);
313  return tmp.advect(time0, time1);
314 }
315 
316 
317 ///////////////////////////////////////////////////////////////////////
318 
319 template<typename GridT, typename InterruptT>
320 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
321  math::TemporalIntegrationScheme TemporalScheme>
322 inline
323 LevelSetMorphing<GridT, InterruptT>::
324 Morph<MapT, SpatialScheme, TemporalScheme>::
325 Morph(LevelSetMorphing<GridT, InterruptT>& parent)
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
337 LevelSetMorphing<GridT, InterruptT>::
338 Morph<MapT, SpatialScheme, TemporalScheme>::
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
352 LevelSetMorphing<GridT, InterruptT>::
353 Morph<MapT, SpatialScheme, TemporalScheme>::
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
451 LevelSetMorphing<GridT, InterruptT>::
452 Morph<MapT, SpatialScheme, TemporalScheme>::
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
481 LevelSetMorphing<GridT, InterruptT>::
482 Morph<MapT, SpatialScheme, TemporalScheme>::
483 sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
484 {
485  using VoxelIterT = typename LeafType::ValueOnCIter;
486  using SamplerT = tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler>;
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
531 LevelSetMorphing<GridT, InterruptT>::
532 Morph<MapT, SpatialScheme, TemporalScheme>::
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
578 LevelSetMorphing<GridT, InterruptT>::
579 Morph<MapT, SpatialScheme, TemporalScheme>::
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
608 LevelSetMorphing<GridT,InterruptT>::
609 Morph<MapT, SpatialScheme, TemporalScheme>::
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;
616  using NumGrad = math::GradientNormSqrd<MapT, SpatialScheme>;
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 
641 ////////////////////////////////////////
642 
643 
644 // Explicit Template Instantiation
645 
646 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
647 
648 #ifdef OPENVDB_INSTANTIATE_LEVELSETMORPH
650 #endif
651 
652 OPENVDB_INSTANTIATE_CLASS LevelSetMorphing<FloatGrid, util::NullInterrupter>;
653 OPENVDB_INSTANTIATE_CLASS LevelSetMorphing<DoubleGrid, util::NullInterrupter>;
654 
655 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
656 
657 
658 } // namespace tools
659 } // namespace OPENVDB_VERSION_NAME
660 } // namespace openvdb
661 
662 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
TemporalIntegrationScheme
Temporal integration schemes.
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:82
GLenum GLint * range
Definition: glcorearb.h:1925
math::Vec3< Real > Vec3R
Definition: Types.h:72
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Return the temporal integration scheme.
void setTemporalScheme(math::TemporalIntegrationScheme s)
Set the spatial finite difference scheme.
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:443
GLboolean invert
Definition: glcorearb.h:549
void setTarget(const GridT &targetGrid)
Redefine the target level set.
Definition: LevelSetMorph.h:74
const GLdouble * v
Definition: glcorearb.h:837
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
GLdouble s
Definition: glad.h:3009
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:239
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
void setAlphaMask(const GridT &maskGrid)
Define the alpha mask.
Definition: LevelSetMorph.h:77
ValueType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
**But if you need a result
Definition: thread.h:613
Performs multi-threaded interface tracking of narrow band level sets.
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:90
int getGrainSize() const
Return the grain size used for multithreading.
math::TemporalIntegrationScheme getTemporalScheme() const
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
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 ...
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
GLdouble n
Definition: glcorearb.h:2008
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition: Math.h:406
void setMaskRange(ValueType min, ValueType max)
Define the range for the (optional) scalar mask.
#define OPENVDB_INSTANTIATE_CLASS
math::BiasedGradientScheme getSpatialScheme() const
Coord Abs(const Coord &xyz)
Definition: Coord.h:517
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=nullptr)
Main constructor.
Definition: LevelSetMorph.h:59
GA_API const UT_StringHolder trans
LevelSetTracker< GridT, InterruptT > TrackerT
Definition: LevelSetMorph.h:52
GLint GLuint mask
Definition: glcorearb.h:124
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
auto get(const UT_ARTIterator< T > &it) -> decltype(it.key())
Definition: UT_ARTMap.h:1073
GLenum target
Definition: glcorearb.h:1667
void setSpatialScheme(math::BiasedGradientScheme s)
Set the spatial finite difference scheme.
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...
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:349
GA_API const UT_StringHolder transform
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:140
GLdouble t
Definition: glad.h:2397
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:85
GLenum mode
Definition: glcorearb.h:99
void setGrainSize(int grainsize)
Set the grain size used for multithreading.
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
static Name mapType()
Return UnitaryMap.
Definition: Maps.h:1717
ValueType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
int getNormCount() const
Return the number of normalizations performed per track or normalize call.
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:87
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:95
size_t advect(ValueType time0, ValueType time1)
Advect the level set from its current time, time0, to its final time, time1. If time0 > time1...
typename LeafManagerType::BufferType BufferType
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:656
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:141
typename LeafManagerType::LeafRange LeafRange
GLboolean r
Definition: glcorearb.h:1222
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition: LevelSetMorph.h:47
GLint GLfloat GLint stencil
Definition: glcorearb.h:1278
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x .
Definition: Math.h:287
__hostdev__ float Sqrt(float x)
Return the square root of a floating-point value.
Definition: NanoVDB.h:1214
GA_API const UT_StringHolder Alpha
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:119
bool ValueType
Definition: NanoVDB.h:5729
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:595
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:337
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:80
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
auto join(It begin, Sentinel end, string_view sep) -> join_view< It, Sentinel >
Definition: format.h:2559
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.