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