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