OpenVDB  3.0.0
LevelSetMorph.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2014 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 ScalarType;
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  ScalarType minMask() const { return mMinMask; }
152 
155  ScalarType maxMask() const { return mDeltaMask + mMinMask; }
156 
164  void setMaskRange(ScalarType min, ScalarType 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(ScalarType time0, ScalarType time1);
183 
184 private:
185 
186  template<math::BiasedGradientScheme SpatialScheme>
187  size_t advect1(ScalarType time0, ScalarType time1);
188 
189  template<math::BiasedGradientScheme SpatialScheme,
190  math::TemporalIntegrationScheme TemporalScheme>
191  size_t advect2(ScalarType time0, ScalarType time1);
192 
193  template<math::BiasedGradientScheme SpatialScheme,
194  math::TemporalIntegrationScheme TemporalScheme,
195  typename MapType>
196  size_t advect3(ScalarType time0, ScalarType time1);
197 
198  TrackerT mTracker;
199  const GridT *mTarget, *mMask;
200  math::BiasedGradientScheme mSpatialScheme;
201  math::TemporalIntegrationScheme mTemporalScheme;
202  ScalarType mMinMask, mDeltaMask;
203  bool mInvertMask;
204 
205  // disallow copy by assignment
206  void operator=(const LevelSetMorphing& other) {}
207 
208  // This templated private class implements all the level set magic.
209  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
210  math::TemporalIntegrationScheme TemporalScheme>
211  class LevelSetMorph
212  {
213  public:
215  LevelSetMorph(LevelSetMorphing<GridT, InterruptT>& parent);
217  LevelSetMorph(const LevelSetMorph& other);
219  LevelSetMorph(LevelSetMorph& other, tbb::split);
221  virtual ~LevelSetMorph() {}
224  size_t advect(ScalarType time0, ScalarType time1);
226  void operator()(const LeafRange& r) const
227  {
228  if (mTask) mTask(const_cast<LevelSetMorph*>(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 LevelSetMorph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
239  private:
240  typedef typename boost::function<void (LevelSetMorph*, const LeafRange&)> FuncType;
241  LevelSetMorphing* mParent;
242  ScalarType mMinAbsS, mMaxAbsS;
243  const MapT* mMap;
244  FuncType mTask;
245 
247  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
248  // method calling tbb
249  void cook(ThreadingMode mode, size_t swapBuffer = 0);
250 
252  typename GridT::ValueType sampleSpeed(ScalarType time0, ScalarType time1, Index speedBuffer);
253  void sampleXformedSpeed(const LeafRange& r, Index speedBuffer);
254  void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer);
255 
256  // Forward Euler advection steps: Phi(result) = Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|;
257  void euler1(const LeafRange& r, ScalarType dt, Index resultBuffer, Index speedBuffer);
258 
259  // Convex combination of Phi and a forward Euler advection steps:
260  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|);
261  void euler2(const LeafRange& r, ScalarType dt, ScalarType alpha,
262  Index phiBuffer, Index resultBuffer, Index speedBuffer);
263 
264  }; // end of private LevelSetMorph class
265 
266 };//end of LevelSetMorphing
267 
268 template<typename GridT, typename InterruptT>
269 inline size_t
270 LevelSetMorphing<GridT, InterruptT>::advect(ScalarType time0, ScalarType time1)
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
292 LevelSetMorphing<GridT, InterruptT>::advect1(ScalarType time0, ScalarType 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  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
311 LevelSetMorphing<GridT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
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
334 LevelSetMorphing<GridT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
335 {
336  LevelSetMorph<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
347 LevelSetMorphing<GridT, InterruptT>::
348 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
349 LevelSetMorph(LevelSetMorphing<GridT, InterruptT>& parent)
350  : mParent(&parent)
351  , mMinAbsS(ScalarType(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
361 LevelSetMorphing<GridT, InterruptT>::
362 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
363 LevelSetMorph(const LevelSetMorph& 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
376 LevelSetMorphing<GridT, InterruptT>::
377 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
378 LevelSetMorph(LevelSetMorph& 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(ScalarType time0, ScalarType 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 ScalarType 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(&LevelSetMorph::euler1, _1, _2, dt, /*result=*/1, /*speed*/2);
411  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
412  this->cook(PARALLEL_FOR, 1);
413  break;
414  case math::TVD_RK2:
415  // Perform one explicit Euler step: t1 = t0 + dt
416  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
417  mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, /*result=*/1, /*speed*/2);
418  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
419  this->cook(PARALLEL_FOR, 1);
420 
421  // Convex combine explict Euler step: t2 = t0 + dt
422  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|)
423  mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.5),
424  /*phi=*/1, /*result=*/1, /*speed*/2);
425  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
426  this->cook(PARALLEL_FOR, 1);
427  break;
428  case math::TVD_RK3:
429  // Perform one explicit Euler step: t1 = t0 + dt
430  // Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]|
431  mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, /*result=*/1, /*speed*/3);
432  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
433  this->cook(PARALLEL_FOR, 1);
434 
435  // Convex combine explict Euler step: t2 = t0 + dt/2
436  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|)
437  mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.75),
438  /*phi=*/1, /*result=*/2, /*speed*/3);
439  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
440  this->cook(PARALLEL_FOR, 2);
441 
442  // Convex combine explict Euler step: t3 = t0 + dt
443  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|)
444  mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(1.0/3.0),
445  /*phi=*/1, /*result=*/2, /*speed*/3);
446  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
447  this->cook(PARALLEL_FOR, 2);
448  break;
449  default:
450  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
451  }//end of compile-time resolved switch
453 
454  time0 += dt;
455  ++countCFL;
456  mParent->mTracker.leafs().removeAuxBuffers();
457 
458  // Track the narrow band
459  mParent->mTracker.track();
460  }//end wile-loop over time
461 
462  return countCFL;//number of CLF propagation steps
463 }
464 
465 template<typename GridT, typename InterruptT>
466 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
467  math::TemporalIntegrationScheme TemporalScheme>
468 inline typename GridT::ValueType
469 LevelSetMorphing<GridT, InterruptT>::
470 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
471 sampleSpeed(ScalarType time0, ScalarType time1, Index speedBuffer)
472 {
473  mMaxAbsS = mMinAbsS;
474  const size_t leafCount = mParent->mTracker.leafs().leafCount();
475  if (leafCount==0 || time0 >= time1) return ScalarType(0);
476 
477  const math::Transform& xform = mParent->mTracker.grid().transform();
478  if (mParent->mTarget->transform() == xform &&
479  (mParent->mMask == NULL || mParent->mMask->transform() == xform)) {
480  mTask = boost::bind(&LevelSetMorph::sampleAlignedSpeed, _1, _2, speedBuffer);
481  } else {
482  mTask = boost::bind(&LevelSetMorph::sampleXformedSpeed, _1, _2, speedBuffer);
483  }
484  this->cook(PARALLEL_REDUCE);
485  if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ScalarType(0);//speed is essentially zero
486  static const ScalarType CFL = (TemporalScheme == math::TVD_RK1 ? ScalarType(0.3) :
487  TemporalScheme == math::TVD_RK2 ? ScalarType(0.9) :
488  ScalarType(1.0))/math::Sqrt(ScalarType(3.0));
489  const ScalarType dt = math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
490  return math::Min(dt, ScalarType(CFL*dx/mMaxAbsS));
491 }
492 
493 template<typename GridT, typename InterruptT>
494 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
495  math::TemporalIntegrationScheme TemporalScheme>
496 inline void
497 LevelSetMorphing<GridT, InterruptT>::
498 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
499 sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
500 {
501  typedef typename LeafType::ValueOnCIter VoxelIterT;
502  typedef tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler> SamplerT;
503  const MapT& map = *mMap;
504  mParent->mTracker.checkInterrupter();
505 
506  typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
507  SamplerT target(targetAcc, mParent->mTarget->transform());
508  if (mParent->mMask == NULL) {
509  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
510  BufferType& speed = leafIter.buffer(speedBuffer);
511  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
512  ScalarType& s = const_cast<ScalarType&>(speed.getValue(voxelIter.pos()));
513  s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
514  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
515  }
516  }
517  } else {
518  const ScalarType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
519  const bool invMask = mParent->isMaskInverted();
520  typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
521  SamplerT mask(maskAcc, mParent->mMask->transform());
522  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
523  BufferType& source = leafIter.buffer(speedBuffer);
524  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
525  const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());//world space
526  const ScalarType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm);
527  ScalarType& s = const_cast<ScalarType&>(source.getValue(voxelIter.pos()));
528  s -= target.wsSample(xyz);
529  s *= invMask ? 1 - a : a;
530  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
531  }
532  }
533  }
534 }
535 
536 template<typename GridT, typename InterruptT>
537 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
538  math::TemporalIntegrationScheme TemporalScheme>
539 inline void
540 LevelSetMorphing<GridT, InterruptT>::
541 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
542 sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
543 {
544  typedef typename LeafType::ValueOnCIter VoxelIterT;
545  mParent->mTracker.checkInterrupter();
546 
547  typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
548 
549  if (mParent->mMask == NULL) {
550  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
551  BufferType& source = leafIter.buffer(speedBuffer);
552  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
553  ScalarType& s = const_cast<ScalarType&>(source.getValue(voxelIter.pos()));
554  s -= target.getValue(voxelIter.getCoord());
555  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
556  }
557  }
558  } else {
559  const ScalarType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
560  const bool invMask = mParent->isMaskInverted();
561  typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
562  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
563  BufferType& source = leafIter.buffer(speedBuffer);
564  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
565  const Coord ijk = voxelIter.getCoord();//index space
566  const ScalarType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm);
567  ScalarType& s = const_cast<ScalarType&>(source.getValue(voxelIter.pos()));
568  s -= target.getValue(ijk);
569  s *= invMask ? 1 - a : a;
570  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
571  }
572  }
573  }
574 }
575 
576 template<typename GridT, typename InterruptT>
577 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
578  math::TemporalIntegrationScheme TemporalScheme>
579 inline void
580 LevelSetMorphing<GridT, InterruptT>::
581 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
582 cook(ThreadingMode mode, size_t swapBuffer)
583 {
584  mParent->mTracker.startInterrupter("Morphing level set");
585 
586  const int grainSize = mParent->mTracker.getGrainSize();
587  const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
588 
589  if (mParent->mTracker.getGrainSize()==0) {
590  (*this)(range);
591  } else if (mode == PARALLEL_FOR) {
592  tbb::parallel_for(range, *this);
593  } else if (mode == PARALLEL_REDUCE) {
594  tbb::parallel_reduce(range, *this);
595  } else {
596  throw std::runtime_error("Undefined threading mode");
597  }
598 
599  mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
600 
601  mParent->mTracker.endInterrupter();
602 }
603 
604 // Forward Euler advection steps:
605 // Phi(result) = Phi(0) - dt * Phi(speed) * |Grad[Phi(0)]|
606 template<typename GridT, typename InterruptT>
607 template <typename MapT,
608  math::BiasedGradientScheme SpatialScheme,
609  math::TemporalIntegrationScheme TemporalScheme>
610 inline void
611 LevelSetMorphing<GridT, InterruptT>::
612 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
613 euler1(const LeafRange& range, ScalarType dt, Index resultBuffer, Index speedBuffer)
614 {
615  typedef math::BIAS_SCHEME<SpatialScheme> SchemeT;
616  typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
617  typedef typename LeafType::ValueOnCIter VoxelIterT;
618  typedef math::GradientNormSqrd<MapT, SpatialScheme> NumGrad;
619 
620  mParent->mTracker.checkInterrupter();
621  const MapT& map = *mMap;
622  StencilT stencil(mParent->mTracker.grid());
623 
624  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
625  BufferType& speed = leafIter.buffer(speedBuffer);
626  BufferType& result = leafIter.buffer(resultBuffer);
627  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
628  const Index n = voxelIter.pos();
629  const ScalarType S = speed.getValue(n);
630  if (math::isApproxZero(S)) continue;
631  stencil.moveTo(voxelIter);
632  result.setValue(n, *voxelIter - dt * S * NumGrad::result(map, stencil));
633  }
634  }
635 }
636 
637 // Convex combination of Phi and a forward Euler advection steps:
638 // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Phi(speed) * |Grad[Phi(0)]|)
639 template<typename GridT, typename InterruptT>
640 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
641  math::TemporalIntegrationScheme TemporalScheme>
642 inline void
643 LevelSetMorphing<GridT, InterruptT>::
644 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
645 euler2(const LeafRange& range, ScalarType dt, ScalarType alpha,
646  Index phiBuffer, Index resultBuffer, Index speedBuffer)
647 {
648  typedef math::BIAS_SCHEME<SpatialScheme> SchemeT;
649  typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
650  typedef typename LeafType::ValueOnCIter VoxelIterT;
651  typedef math::GradientNormSqrd<MapT, SpatialScheme> NumGrad;
652 
653  mParent->mTracker.checkInterrupter();
654  const MapT& map = *mMap;
655  const ScalarType beta = ScalarType(1.0) - alpha;
656  StencilT stencil(mParent->mTracker.grid());
657 
658  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
659  BufferType& speed = leafIter.buffer(speedBuffer);
660  BufferType& result = leafIter.buffer(resultBuffer);
661  BufferType& phi = leafIter.buffer(phiBuffer);
662  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
663  const Index n = voxelIter.pos();
664  const ScalarType S = speed.getValue(n);
665  if (math::isApproxZero(S)) continue;
666  stencil.moveTo(voxelIter);
667  const ScalarType G = NumGrad::result(map, stencil);
668  result.setValue(n, alpha * phi.getValue(n) + beta * (*voxelIter - dt * S * G));
669  }
670  }
671 }
672 
673 } // namespace tools
674 } // namespace OPENVDB_VERSION_NAME
675 } // namespace openvdb
676 
677 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
678 
679 // Copyright (c) 2012-2014 DreamWorks Animation LLC
680 // All rights reserved. This software is distributed under the
681 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:111
Definition: FiniteDifference.h:265
void setMaskRange(ScalarType min, ScalarType max)
Define the range for the (optional) scalar mask.
Definition: LevelSetMorph.h:164
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:284
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:119
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
TrackerT::BufferType BufferType
Definition: LevelSetMorph.h:82
Name mapType() const
Return the transformation map's type-name.
Definition: Transform.h:93
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:116
TrackerT::LeafType LeafType
Definition: LevelSetMorph.h:81
void setGrainSize(int grainsize)
Set the grain size used for multithreading.
Definition: LevelSetMorph.h:147
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:129
int getGrainSize() const
Return the grain size used for multithreading.
Definition: LevelSetMorph.h:144
TrackerT::ValueType ScalarType
Definition: LevelSetMorph.h:83
ScalarType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetMorph.h:151
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
size_t advect(ScalarType time0, ScalarType time1)
Advect the level set from its current time, time0, to its final time, time1. If time0 > time1...
Definition: LevelSetMorph.h:270
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition: LevelSetMorph.h:74
TreeType::LeafNodeType LeafType
Definition: LevelSetTracker.h:73
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:606
GridT::TreeType TreeType
Definition: LevelSetMorph.h:78
LevelSetTracker< GridT, InterruptT > TrackerT
Definition: LevelSetMorph.h:79
void setAlphaMask(const GridT &maskGrid)
Define the alpha mask.
Definition: LevelSetMorph.h:106
void setTarget(const GridT &targetGrid)
Redefine the target level set.
Definition: LevelSetMorph.h:103
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:693
TrackerT::LeafRange LeafRange
Definition: LevelSetMorph.h:80
#define OPENVDB_VERSION_NAME
Definition: version.h:43
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:141
Index32 Index
Definition: Types.h:59
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
math::Vec3< Real > Vec3R
Definition: Types.h:77
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:263
Definition: Exceptions.h:39
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:121
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:114
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:326
Definition: FiniteDifference.h:264
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:545
Definition: Exceptions.h:88
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:122
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=NULL)
Main constructor.
Definition: LevelSetMorph.h:86
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
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
ScalarType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetMorph.h:155
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:68
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:124
Definition: FiniteDifference.h:198
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
LeafManagerType::BufferType BufferType
Definition: LevelSetTracker.h:78
Definition: FiniteDifference.h:194
GridT GridType
Definition: LevelSetMorph.h:77
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:314
Definition: FiniteDifference.h:263
TreeType::ValueType ValueType
Definition: LevelSetTracker.h:74
int getNormCount() const
Return the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:139
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:109
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:134
virtual ~LevelSetMorphing()
Definition: LevelSetMorph.h:100
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261