HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LevelSetMorph.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2018 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 //
29 ///////////////////////////////////////////////////////////////////////////
30 
31 /// @author Ken Museth
32 ///
33 /// @file tools/LevelSetMorph.h
34 ///
35 /// @brief Shape morphology of level sets. Morphing from a source
36 /// narrow-band level sets to a target narrow-band level set.
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 #include <functional>
45 #include <limits>
46 
47 
48 namespace openvdb {
50 namespace OPENVDB_VERSION_NAME {
51 namespace tools {
52 
53 /// @brief Shape morphology of level sets. Morphing from a source
54 /// narrow-band level sets to a target narrow-band level set.
55 ///
56 /// @details
57 /// The @c InterruptType template argument below refers to any class
58 /// with the following interface:
59 /// @code
60 /// class Interrupter {
61 /// ...
62 /// public:
63 /// void start(const char* name = nullptr) // called when computations begin
64 /// void end() // called when computations end
65 /// bool wasInterrupted(int percent=-1) // return true to break computation
66 /// };
67 /// @endcode
68 ///
69 /// @note If no template argument is provided for this InterruptType,
70 /// the util::NullInterrupter is used, which implies that all interrupter
71 /// calls are no-ops (i.e., they incur no computational overhead).
72 template<typename GridT,
73  typename InterruptT = util::NullInterrupter>
75 {
76 public:
77  using GridType = GridT;
78  using TreeType = typename GridT::TreeType;
80  using LeafRange = typename TrackerT::LeafRange;
81  using LeafType = typename TrackerT::LeafType;
82  using BufferType = typename TrackerT::BufferType;
83  using ValueType = typename TrackerT::ValueType;
84 
85  /// Main constructor
86  LevelSetMorphing(GridT& sourceGrid, const GridT& targetGrid, InterruptT* interrupt = nullptr)
87  : mTracker(sourceGrid, interrupt)
88  , mTarget(&targetGrid)
89  , mMask(nullptr)
90  , mSpatialScheme(math::HJWENO5_BIAS)
91  , mTemporalScheme(math::TVD_RK2)
92  , mMinMask(0)
93  , mDeltaMask(1)
94  , mInvertMask(false)
95  {
96  }
97 
98  virtual ~LevelSetMorphing() {}
99 
100  /// Redefine the target level set
101  void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; }
102 
103  /// Define the alpha mask
104  void setAlphaMask(const GridT& maskGrid) { mMask = &maskGrid; }
105 
106  /// Return the spatial finite-difference scheme
107  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
108  /// Set the spatial finite-difference scheme
109  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
110 
111  /// Return the temporal integration scheme
112  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
113  /// Set the temporal integration scheme
114  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
115 
116  /// Return the spatial finite-difference scheme
118  {
119  return mTracker.getSpatialScheme();
120  }
121  /// Set the spatial finite-difference scheme
123  {
124  mTracker.setSpatialScheme(scheme);
125  }
126  /// Return the temporal integration scheme
128  {
129  return mTracker.getTemporalScheme();
130  }
131  /// Set the temporal integration scheme
133  {
134  mTracker.setTemporalScheme(scheme);
135  }
136  /// Return the number of normalizations performed per track or normalize call.
137  int getNormCount() const { return mTracker.getNormCount(); }
138  /// Set the number of normalizations performed per track or normalize call.
139  void setNormCount(int n) { mTracker.setNormCount(n); }
140 
141  /// Return the grain size used for multithreading
142  int getGrainSize() const { return mTracker.getGrainSize(); }
143  /// @brief Set the grain size used for multithreading.
144  /// @note A grain size of 0 or less disables multithreading!
145  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
146 
147  /// @brief Return the minimum value of the mask to be used for the
148  /// derivation of a smooth alpha value.
149  ValueType minMask() const { return mMinMask; }
150 
151  /// @brief Return the maximum value of the mask to be used for the
152  /// derivation of a smooth alpha value.
153  ValueType maxMask() const { return mDeltaMask + mMinMask; }
154 
155  /// @brief Define the range for the (optional) scalar mask.
156  /// @param min Minimum value of the range.
157  /// @param max Maximum value of the range.
158  /// @details Mask values outside the range maps to alpha values of
159  /// respectfully zero and one, and values inside the range maps
160  /// smoothly to 0->1 (unless of course the mask is inverted).
161  /// @throw ValueError if @a min is not smaller than @a max.
163  {
164  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
165  mMinMask = min;
166  mDeltaMask = max-min;
167  }
168 
169  /// @brief Return true if the mask is inverted, i.e. min->max in the
170  /// original mask maps to 1->0 in the inverted alpha mask.
171  bool isMaskInverted() const { return mInvertMask; }
172  /// @brief Invert the optional mask, i.e. min->max in the original
173  /// mask maps to 1->0 in the inverted alpha mask.
174  void invertMask(bool invert=true) { mInvertMask = invert; }
175 
176  /// @brief Advect the level set from its current time, @a time0, to its
177  /// final time, @a time1. If @a time0 > @a time1, perform backward advection.
178  ///
179  /// @return the number of CFL iterations used to advect from @a time0 to @a time1
180  size_t advect(ValueType time0, ValueType time1);
181 
182 private:
183 
184  // disallow copy construction and copy by assignment!
185  LevelSetMorphing(const LevelSetMorphing&);// not implemented
186  LevelSetMorphing& operator=(const LevelSetMorphing&);// not implemented
187 
188  template<math::BiasedGradientScheme SpatialScheme>
189  size_t advect1(ValueType time0, ValueType time1);
190 
191  template<math::BiasedGradientScheme SpatialScheme,
192  math::TemporalIntegrationScheme TemporalScheme>
193  size_t advect2(ValueType time0, ValueType time1);
194 
195  template<math::BiasedGradientScheme SpatialScheme,
196  math::TemporalIntegrationScheme TemporalScheme,
197  typename MapType>
198  size_t advect3(ValueType time0, ValueType time1);
199 
200  TrackerT mTracker;
201  const GridT *mTarget, *mMask;
202  math::BiasedGradientScheme mSpatialScheme;
203  math::TemporalIntegrationScheme mTemporalScheme;
204  ValueType mMinMask, mDeltaMask;
205  bool mInvertMask;
206 
207  // This templated private class implements all the level set magic.
208  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
209  math::TemporalIntegrationScheme TemporalScheme>
210  struct Morph
211  {
212  /// Main constructor
214  /// Shallow copy constructor called by tbb::parallel_for() threads
215  Morph(const Morph& other);
216  /// Shallow copy constructor called by tbb::parallel_reduce() threads
217  Morph(Morph& other, tbb::split);
218  /// destructor
219  virtual ~Morph() {}
220  /// Advect the level set from its current time, time0, to its final time, time1.
221  /// @return number of CFL iterations
222  size_t advect(ValueType time0, ValueType time1);
223  /// Used internally by tbb::parallel_for()
224  void operator()(const LeafRange& r) const
225  {
226  if (mTask) mTask(const_cast<Morph*>(this), r);
227  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
228  }
229  /// Used internally by tbb::parallel_reduce()
230  void operator()(const LeafRange& r)
231  {
232  if (mTask) mTask(this, r);
233  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
234  }
235  /// This is only called by tbb::parallel_reduce() threads
236  void join(const Morph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
237 
238  /// Enum to define the type of multithreading
239  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
240  // method calling tbb
241  void cook(ThreadingMode mode, size_t swapBuffer = 0);
242 
243  /// Sample field and return the CFT time step
244  typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer);
245  void sampleXformedSpeed(const LeafRange& r, Index speedBuffer);
246  void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer);
247 
248  // Convex combination of Phi and a forward Euler advection steps:
249  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|);
250  template <int Nominator, int Denominator>
251  void euler(const LeafRange&, ValueType, Index, Index, Index);
252  inline void euler01(const LeafRange& r, ValueType t, Index s) {this->euler<0,1>(r,t,0,1,s);}
253  inline void euler12(const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
254  inline void euler34(const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
255  inline void euler13(const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
256 
257  using FuncType = typename std::function<void (Morph*, const LeafRange&)>;
258  LevelSetMorphing* mParent;
259  ValueType mMinAbsS, mMaxAbsS;
260  const MapT* mMap;
261  FuncType mTask;
262  }; // end of private Morph struct
263 
264 };//end of LevelSetMorphing
265 
266 template<typename GridT, typename InterruptT>
267 inline size_t
269 {
270  switch (mSpatialScheme) {
271  case math::FIRST_BIAS:
272  return this->advect1<math::FIRST_BIAS >(time0, time1);
273  //case math::SECOND_BIAS:
274  //return this->advect1<math::SECOND_BIAS >(time0, time1);
275  //case math::THIRD_BIAS:
276  //return this->advect1<math::THIRD_BIAS >(time0, time1);
277  //case math::WENO5_BIAS:
278  //return this->advect1<math::WENO5_BIAS >(time0, time1);
279  case math::HJWENO5_BIAS:
280  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
281  case math::SECOND_BIAS:
282  case math::THIRD_BIAS:
283  case math::WENO5_BIAS:
284  case math::UNKNOWN_BIAS:
285  default:
286  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
287  }
288  return 0;
289 }
290 
291 template<typename GridT, typename InterruptT>
292 template<math::BiasedGradientScheme SpatialScheme>
293 inline size_t
294 LevelSetMorphing<GridT, InterruptT>::advect1(ValueType time0, ValueType time1)
295 {
296  switch (mTemporalScheme) {
297  case math::TVD_RK1:
298  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
299  case math::TVD_RK2:
300  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
301  case math::TVD_RK3:
302  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
303  case math::UNKNOWN_TIS:
304  default:
305  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
306  }
307  return 0;
308 }
309 
310 template<typename GridT, typename InterruptT>
311 template<math::BiasedGradientScheme SpatialScheme,
312  math::TemporalIntegrationScheme TemporalScheme>
313 inline size_t
314 LevelSetMorphing<GridT, InterruptT>::advect2(ValueType time0, ValueType time1)
315 {
316  const math::Transform& trans = mTracker.grid().transform();
317  if (trans.mapType() == math::UniformScaleMap::mapType()) {
318  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
319  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
320  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
321  time0, time1);
322  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
323  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
324  } else if (trans.mapType() == math::TranslationMap::mapType()) {
325  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
326  } else {
327  OPENVDB_THROW(ValueError, "MapType not supported!");
328  }
329  return 0;
330 }
331 
332 template<typename GridT, typename InterruptT>
333 template<math::BiasedGradientScheme SpatialScheme,
334  math::TemporalIntegrationScheme TemporalScheme,
335  typename MapT>
336 inline size_t
337 LevelSetMorphing<GridT, InterruptT>::advect3(ValueType time0, ValueType time1)
338 {
339  Morph<MapT, SpatialScheme, TemporalScheme> tmp(*this);
340  return tmp.advect(time0, time1);
341 }
342 
343 
344 ///////////////////////////////////////////////////////////////////////
345 
346 template<typename GridT, typename InterruptT>
347 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
348  math::TemporalIntegrationScheme TemporalScheme>
349 inline
350 LevelSetMorphing<GridT, InterruptT>::
351 Morph<MapT, SpatialScheme, TemporalScheme>::
352 Morph(LevelSetMorphing<GridT, InterruptT>& parent)
353  : mParent(&parent)
354  , mMinAbsS(ValueType(1e-6))
355  , mMap(parent.mTracker.grid().transform().template constMap<MapT>().get())
356  , mTask(0)
357 {
358 }
359 
360 template<typename GridT, typename InterruptT>
361 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
362  math::TemporalIntegrationScheme TemporalScheme>
363 inline
364 LevelSetMorphing<GridT, InterruptT>::
365 Morph<MapT, SpatialScheme, TemporalScheme>::
366 Morph(const Morph& other)
367  : mParent(other.mParent)
368  , mMinAbsS(other.mMinAbsS)
369  , mMaxAbsS(other.mMaxAbsS)
370  , mMap(other.mMap)
371  , mTask(other.mTask)
372 {
373 }
374 
375 template<typename GridT, typename InterruptT>
376 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
377  math::TemporalIntegrationScheme TemporalScheme>
378 inline
379 LevelSetMorphing<GridT, InterruptT>::
380 Morph<MapT, SpatialScheme, TemporalScheme>::
381 Morph(Morph& other, tbb::split)
382  : mParent(other.mParent)
383  , mMinAbsS(other.mMinAbsS)
384  , mMaxAbsS(other.mMaxAbsS)
385  , mMap(other.mMap)
386  , mTask(other.mTask)
387 {
388 }
389 
390 template<typename GridT, typename InterruptT>
391 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
392  math::TemporalIntegrationScheme TemporalScheme>
393 inline size_t
396 advect(ValueType time0, ValueType time1)
397 {
398  namespace ph = std::placeholders;
399 
400  // Make sure we have enough temporal auxiliary buffers for the time
401  // integration AS WELL AS an extra buffer with the speed function!
402  static const Index auxBuffers = 1 + (TemporalScheme == math::TVD_RK3 ? 2 : 1);
403  size_t countCFL = 0;
404  while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
405  mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
406 
407  const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
408  if ( math::isZero(dt) ) break;//V is essentially zero so terminate
409 
410  OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time
411  switch(TemporalScheme) {
412  case math::TVD_RK1:
413  // Perform one explicit Euler step: t1 = t0 + dt
414  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
415  mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2);
416 
417  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
418  this->cook(PARALLEL_FOR, 1);
419  break;
420  case math::TVD_RK2:
421  // Perform one explicit Euler step: t1 = t0 + dt
422  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
423  mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2);
424 
425  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
426  this->cook(PARALLEL_FOR, 1);
427 
428  // Convex combine explict Euler step: t2 = t0 + dt
429  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|)
430  mTask = std::bind(&Morph::euler12, ph::_1, ph::_2, dt);
431 
432  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
433  this->cook(PARALLEL_FOR, 1);
434  break;
435  case math::TVD_RK3:
436  // Perform one explicit Euler step: t1 = t0 + dt
437  // Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]|
438  mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/3);
439 
440  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
441  this->cook(PARALLEL_FOR, 1);
442 
443  // Convex combine explict Euler step: t2 = t0 + dt/2
444  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|)
445  mTask = std::bind(&Morph::euler34, ph::_1, ph::_2, dt);
446 
447  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
448  this->cook(PARALLEL_FOR, 2);
449 
450  // Convex combine explict Euler step: t3 = t0 + dt
451  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|)
452  mTask = std::bind(&Morph::euler13, ph::_1, ph::_2, dt);
453 
454  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
455  this->cook(PARALLEL_FOR, 2);
456  break;
457  case math::UNKNOWN_TIS:
458  default:
459  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
460  }//end of compile-time resolved switch
462 
463  time0 += dt;
464  ++countCFL;
465  mParent->mTracker.leafs().removeAuxBuffers();
466 
467  // Track the narrow band
468  mParent->mTracker.track();
469  }//end wile-loop over time
470 
471  return countCFL;//number of CLF propagation steps
472 }
473 
474 template<typename GridT, typename InterruptT>
475 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
476  math::TemporalIntegrationScheme TemporalScheme>
477 inline typename GridT::ValueType
478 LevelSetMorphing<GridT, InterruptT>::
479 Morph<MapT, SpatialScheme, TemporalScheme>::
480 sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer)
481 {
482  namespace ph = std::placeholders;
483 
484  mMaxAbsS = mMinAbsS;
485  const size_t leafCount = mParent->mTracker.leafs().leafCount();
486  if (leafCount==0 || time0 >= time1) return ValueType(0);
487 
488  const math::Transform& xform = mParent->mTracker.grid().transform();
489  if (mParent->mTarget->transform() == xform &&
490  (mParent->mMask == nullptr || mParent->mMask->transform() == xform)) {
491  mTask = std::bind(&Morph::sampleAlignedSpeed, ph::_1, ph::_2, speedBuffer);
492  } else {
493  mTask = std::bind(&Morph::sampleXformedSpeed, ph::_1, ph::_2, speedBuffer);
494  }
495  this->cook(PARALLEL_REDUCE);
496  if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ValueType(0);//speed is essentially zero
497  static const ValueType CFL = (TemporalScheme == math::TVD_RK1 ? ValueType(0.3) :
498  TemporalScheme == math::TVD_RK2 ? ValueType(0.9) :
499  ValueType(1.0))/math::Sqrt(ValueType(3.0));
500  const ValueType dt = math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
501  return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
502 }
503 
504 template<typename GridT, typename InterruptT>
505 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
506  math::TemporalIntegrationScheme TemporalScheme>
507 inline void
508 LevelSetMorphing<GridT, InterruptT>::
509 Morph<MapT, SpatialScheme, TemporalScheme>::
510 sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
511 {
512  using VoxelIterT = typename LeafType::ValueOnCIter;
513  using SamplerT = tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler>;
514 
515  const MapT& map = *mMap;
516  mParent->mTracker.checkInterrupter();
517 
518  typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
519  SamplerT target(targetAcc, mParent->mTarget->transform());
520  if (mParent->mMask == nullptr) {
521  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
522  ValueType* speed = leafIter.buffer(speedBuffer).data();
523  bool isZero = true;
524  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
525  ValueType& s = speed[voxelIter.pos()];
526  s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
527  if (!math::isApproxZero(s)) isZero = false;
528  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
529  }
530  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
531  }
532  } else {
533  const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
534  const bool invMask = mParent->isMaskInverted();
535  typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
536  SamplerT mask(maskAcc, mParent->mMask->transform());
537  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
538  ValueType* speed = leafIter.buffer(speedBuffer).data();
539  bool isZero = true;
540  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
541  const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());//world space
542  const ValueType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm);
543  ValueType& s = speed[voxelIter.pos()];
544  s -= target.wsSample(xyz);
545  s *= invMask ? 1 - a : a;
546  if (!math::isApproxZero(s)) isZero = false;
547  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
548  }
549  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
550  }
551  }
552 }
553 
554 template<typename GridT, typename InterruptT>
555 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
556  math::TemporalIntegrationScheme TemporalScheme>
557 inline void
558 LevelSetMorphing<GridT, InterruptT>::
559 Morph<MapT, SpatialScheme, TemporalScheme>::
560 sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
561 {
562  using VoxelIterT = typename LeafType::ValueOnCIter;
563 
564  mParent->mTracker.checkInterrupter();
565 
566  typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
567 
568  if (mParent->mMask == nullptr) {
569  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
570  ValueType* speed = leafIter.buffer(speedBuffer).data();
571  bool isZero = true;
572  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
573  ValueType& s = speed[voxelIter.pos()];
574  s -= target.getValue(voxelIter.getCoord());
575  if (!math::isApproxZero(s)) isZero = false;
576  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
577  }
578  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
579  }
580  } else {
581  const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
582  const bool invMask = mParent->isMaskInverted();
583  typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
584  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
585  ValueType* speed = leafIter.buffer(speedBuffer).data();
586  bool isZero = true;
587  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
588  const Coord ijk = voxelIter.getCoord();//index space
589  const ValueType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm);
590  ValueType& s = speed[voxelIter.pos()];
591  s -= target.getValue(ijk);
592  s *= invMask ? 1 - a : a;
593  if (!math::isApproxZero(s)) isZero = false;
594  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
595  }
596  if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
597  }
598  }
599 }
600 
601 template<typename GridT, typename InterruptT>
602 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
603  math::TemporalIntegrationScheme TemporalScheme>
604 inline void
605 LevelSetMorphing<GridT, InterruptT>::
606 Morph<MapT, SpatialScheme, TemporalScheme>::
607 cook(ThreadingMode mode, size_t swapBuffer)
608 {
609  mParent->mTracker.startInterrupter("Morphing level set");
610 
611  const int grainSize = mParent->mTracker.getGrainSize();
612  const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
613 
614  if (mParent->mTracker.getGrainSize()==0) {
615  (*this)(range);
616  } else if (mode == PARALLEL_FOR) {
617  tbb::parallel_for(range, *this);
618  } else if (mode == PARALLEL_REDUCE) {
619  tbb::parallel_reduce(range, *this);
620  } else {
621  OPENVDB_THROW(ValueError, "expected threading mode " << int(PARALLEL_FOR)
622  << " or " << int(PARALLEL_REDUCE) << ", got " << int(mode));
623  }
624 
625  mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
626 
627  mParent->mTracker.endInterrupter();
628 }
629 
630 template<typename GridT, typename InterruptT>
631 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
632  math::TemporalIntegrationScheme TemporalScheme>
633 template <int Nominator, int Denominator>
634 inline void
635 LevelSetMorphing<GridT,InterruptT>::
636 Morph<MapT, SpatialScheme, TemporalScheme>::
637 euler(const LeafRange& range, ValueType dt,
638  Index phiBuffer, Index resultBuffer, Index speedBuffer)
639 {
640  using SchemeT = math::BIAS_SCHEME<SpatialScheme>;
641  using StencilT = typename SchemeT::template ISStencil<GridType>::StencilType;
642  using VoxelIterT = typename LeafType::ValueOnCIter;
643  using NumGrad = math::GradientNormSqrd<MapT, SpatialScheme>;
644 
645  static const ValueType Alpha = ValueType(Nominator)/ValueType(Denominator);
646  static const ValueType Beta = ValueType(1) - Alpha;
647 
648  mParent->mTracker.checkInterrupter();
649  const MapT& map = *mMap;
650  StencilT stencil(mParent->mTracker.grid());
651 
652  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
653  const ValueType* speed = leafIter.buffer(speedBuffer).data();
655  const ValueType* phi = leafIter.buffer(phiBuffer).data();
656  ValueType* result = leafIter.buffer(resultBuffer).data();
657  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
658  const Index n = voxelIter.pos();
659  if (math::isApproxZero(speed[n])) continue;
660  stencil.moveTo(voxelIter);
661  const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
662  result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
663  }//loop over active voxels in the leaf of the mask
664  }//loop over leafs of the level set
665 }
666 
667 } // namespace tools
668 } // namespace OPENVDB_VERSION_NAME
669 } // namespace openvdb
670 
671 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
672 
673 // Copyright (c) 2012-2018 DreamWorks Animation LLC
674 // All rights reserved. This software is distributed under the
675 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
TemporalIntegrationScheme
Temporal integration schemes.
GLenum GLint * range
Definition: glcorearb.h:1924
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
math::Vec3< Real > Vec3R
Definition: Types.h:79
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:395
GLboolean invert
Definition: glcorearb.h:548
void setTarget(const GridT &targetGrid)
Redefine the target level set.
const GLdouble * v
Definition: glcorearb.h:836
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
GLint GLuint mask
Definition: glcorearb.h:123
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
void setAlphaMask(const GridT &maskGrid)
Define the alpha mask.
ValueType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Performs multi-threaded interface tracking of narrow band level sets.
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
int getGrainSize() const
Return the grain size used for multithreading.
math::TemporalIntegrationScheme getTemporalScheme() const
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
png_infop png_bytep * trans
Definition: png.h:2520
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:2007
void setMaskRange(ValueType min, ValueType max)
Define the range for the (optional) scalar mask.
math::BiasedGradientScheme getSpatialScheme() const
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=nullptr)
Main constructor.
Definition: LevelSetMorph.h:86
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
LevelSetTracker< GridT, InterruptT > TrackerT
Definition: LevelSetMorph.h:79
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:715
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
GLenum target
Definition: glcorearb.h:1666
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:320
GA_API const UT_StringHolder transform
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:129
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:358
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
GLint GLfloat GLint stencil
Definition: glcorearb.h:1277
GLenum mode
Definition: glcorearb.h:98
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:1731
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.
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
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
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:610
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:130
typename LeafManagerType::LeafRange LeafRange
GLboolean r
Definition: glcorearb.h:1221
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition: LevelSetMorph.h:74
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x .
Definition: Math.h:256
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
GA_API const UT_StringHolder Alpha
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:549
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:308
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.