OpenVDB  2.3.0
LevelSetMorph.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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, const GridT& targetGrid, InterruptT* interrupt = NULL):
87  mTracker(sourceGrid, interrupt), mTarget(&targetGrid),
88  mSpatialScheme(math::HJWENO5_BIAS),
89  mTemporalScheme(math::TVD_RK2) {}
90 
91  virtual ~LevelSetMorphing() {};
92 
94  void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; }
95 
97  math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
99  void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
100 
102  math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
104  void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
105 
108  {
109  return mTracker.getSpatialScheme();
110  }
113  {
114  mTracker.setSpatialScheme(scheme);
115  }
118  {
119  return mTracker.getTemporalScheme();
120  }
123  {
124  mTracker.setTemporalScheme(scheme);
125  }
127  int getNormCount() const { return mTracker.getNormCount(); }
129  void setNormCount(int n) { mTracker.setNormCount(n); }
130 
132  int getGrainSize() const { return mTracker.getGrainSize(); }
135  void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
136 
141  size_t advect(ScalarType time0, ScalarType time1);
142 
143 private:
144 
145  // This templated private class implements all the level set magic.
146  template<typename MapT, math::BiasedGradientScheme SpatialScheme,
147  math::TemporalIntegrationScheme TemporalScheme>
148  class LevelSetMorph
149  {
150  public:
152  LevelSetMorph(TrackerT& tracker, const GridT* target);
154  LevelSetMorph(const LevelSetMorph& other);
156  LevelSetMorph(LevelSetMorph& other, tbb::split);
158  virtual ~LevelSetMorph() {}
161  size_t advect(ScalarType time0, ScalarType time1);
163  void operator()(const LeafRange& r) const
164  {
165  if (mTask) mTask(const_cast<LevelSetMorph*>(this), r);
166  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
167  }
169  void operator()(const LeafRange& r)
170  {
171  if (mTask) mTask(this, r);
172  else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
173  }
175  void join(const LevelSetMorph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
176  private:
177  typedef typename boost::function<void (LevelSetMorph*, const LeafRange&)> FuncType;
178  TrackerT* mTracker;
179  const GridT* mTarget;
180  ScalarType mMinAbsS, mMaxAbsS;
181  const MapT* mMap;
182  FuncType mTask;
183 
185  enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
186  // method calling tbb
187  void cook(ThreadingMode mode, size_t swapBuffer = 0);
188 
190  typename GridT::ValueType sampleSpeed(
191  ScalarType time0, ScalarType time1, Index speedBuffer);
192  void sampleXformedSpeed(const LeafRange& r, Index speedBuffer);
193  void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer);
194 
195  // Forward Euler advection steps: Phi(result) = Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|;
196  void euler1(const LeafRange& r, ScalarType dt, Index resultBuffer, Index speedBuffer);
197 
198  // Convex combination of Phi and a forward Euler advection steps:
199  // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|);
200  void euler2(const LeafRange& r, ScalarType dt, ScalarType alpha,
201  Index phiBuffer, Index resultBuffer, Index speedBuffer);
202 
203  }; // end of private LevelSetMorph class
204 
205  template<math::BiasedGradientScheme SpatialScheme>
206  size_t advect1(ScalarType time0, ScalarType time1);
207 
208  template<math::BiasedGradientScheme SpatialScheme,
209  math::TemporalIntegrationScheme TemporalScheme>
210  size_t advect2(ScalarType time0, ScalarType time1);
211 
212  template<math::BiasedGradientScheme SpatialScheme,
213  math::TemporalIntegrationScheme TemporalScheme,
214  typename MapType>
215  size_t advect3(ScalarType time0, ScalarType time1);
216 
217  TrackerT mTracker;
218  const GridT* mTarget;
219  math::BiasedGradientScheme mSpatialScheme;
220  math::TemporalIntegrationScheme mTemporalScheme;
221 
222  // disallow copy by assignment
223  void operator=(const LevelSetMorphing& other) {}
224 
225 };//end of LevelSetMorphing
226 
227 template<typename GridT, typename InterruptT>
228 inline size_t
230 {
231  switch (mSpatialScheme) {
232  case math::FIRST_BIAS:
233  return this->advect1<math::FIRST_BIAS >(time0, time1);
234  //case math::SECOND_BIAS:
235  //return this->advect1<math::SECOND_BIAS >(time0, time1);
236  //case math::THIRD_BIAS:
237  //return this->advect1<math::THIRD_BIAS >(time0, time1);
238  //case math::WENO5_BIAS:
239  //return this->advect1<math::WENO5_BIAS >(time0, time1);
240  case math::HJWENO5_BIAS:
241  return this->advect1<math::HJWENO5_BIAS>(time0, time1);
242  default:
243  OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
244  }
245  return 0;
246 }
247 
248 template<typename GridT, typename InterruptT>
249 template<math::BiasedGradientScheme SpatialScheme>
250 inline size_t
251 LevelSetMorphing<GridT, InterruptT>::advect1(ScalarType time0, ScalarType time1)
252 {
253  switch (mTemporalScheme) {
254  case math::TVD_RK1:
255  return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
256  case math::TVD_RK2:
257  return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
258  case math::TVD_RK3:
259  return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
260  default:
261  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
262  }
263  return 0;
264 }
265 
266 template<typename GridT, typename InterruptT>
267 template<math::BiasedGradientScheme SpatialScheme,
268  math::TemporalIntegrationScheme TemporalScheme>
269 inline size_t
270 LevelSetMorphing<GridT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
271 {
272  const math::Transform& trans = mTracker.grid().transform();
273  if (trans.mapType() == math::UniformScaleMap::mapType()) {
274  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
275  } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
276  return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
277  time0, time1);
278  } else if (trans.mapType() == math::UnitaryMap::mapType()) {
279  return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
280  } else if (trans.mapType() == math::TranslationMap::mapType()) {
281  return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
282  } else {
283  OPENVDB_THROW(ValueError, "MapType not supported!");
284  }
285  return 0;
286 }
287 
288 template<typename GridT, typename InterruptT>
289 template<math::BiasedGradientScheme SpatialScheme,
290  math::TemporalIntegrationScheme TemporalScheme,
291  typename MapT>
292 inline size_t
293 LevelSetMorphing<GridT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
294 {
295  LevelSetMorph<MapT, SpatialScheme, TemporalScheme> tmp(mTracker, mTarget);
296  return tmp.advect(time0, time1);
297 }
298 
299 
301 
302 
303 template<typename GridT, typename InterruptT>
304 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
305  math::TemporalIntegrationScheme TemporalScheme>
306 inline
307 LevelSetMorphing<GridT, InterruptT>::
308 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
309 LevelSetMorph(TrackerT& tracker, const GridT* target):
310  mTracker(&tracker),
311  mTarget(target),
312  mMinAbsS(1e-6),
313  mMap(tracker.grid().transform().template constMap<MapT>().get()),
314  mTask(0)
315 {
316 }
317 
318 template<typename GridT, typename InterruptT>
319 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
320  math::TemporalIntegrationScheme TemporalScheme>
321 inline
322 LevelSetMorphing<GridT, InterruptT>::
323 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
324 LevelSetMorph(const LevelSetMorph& other):
325  mTracker(other.mTracker),
326  mTarget(other.mTarget),
327  mMinAbsS(other.mMinAbsS),
328  mMaxAbsS(other.mMaxAbsS),
329  mMap(other.mMap),
330  mTask(other.mTask)
331 {
332 }
333 
334 template<typename GridT, typename InterruptT>
335 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
336  math::TemporalIntegrationScheme TemporalScheme>
337 inline
338 LevelSetMorphing<GridT, InterruptT>::
339 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
340 LevelSetMorph(LevelSetMorph& other, tbb::split):
341  mTracker(other.mTracker),
342  mTarget(other.mTarget),
343  mMinAbsS(other.mMinAbsS),
344  mMaxAbsS(other.mMaxAbsS),
345  mMap(other.mMap),
346  mTask(other.mTask)
347 {
348 }
349 
350 template<typename GridT, typename InterruptT>
351 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
352  math::TemporalIntegrationScheme TemporalScheme>
353 inline size_t
356 advect(ScalarType time0, ScalarType time1)
357 {
358  // Make sure we have enough temporal auxiliary buffers for the time
359  // integration AS WELL AS an extra buffer with the speed function!
360  static const Index auxBuffers = 1 + (TemporalScheme == math::TVD_RK3 ? 2 : 1);
361  size_t countCFL = 0;
362  while (time0 < time1 && mTracker->checkInterrupter()) {
363  mTracker->leafs().rebuildAuxBuffers(auxBuffers);
364 
365  const ScalarType dt = this->sampleSpeed(time0, time1, auxBuffers);
366  if ( math::isZero(dt) ) break;//V is essentially zero so terminate
367 
368  OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time
369  switch(TemporalScheme) {
370  case math::TVD_RK1:
371  // Perform one explicit Euler step: t1 = t0 + dt
372  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
373  mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, /*result=*/1, /*speed*/2);
374  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
375  this->cook(PARALLEL_FOR, 1);
376  break;
377  case math::TVD_RK2:
378  // Perform one explicit Euler step: t1 = t0 + dt
379  // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
380  mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, /*result=*/1, /*speed*/2);
381  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
382  this->cook(PARALLEL_FOR, 1);
383 
384  // Convex combine explict Euler step: t2 = t0 + dt
385  // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|)
386  mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.5),
387  /*phi=*/1, /*result=*/1, /*speed*/2);
388  // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
389  this->cook(PARALLEL_FOR, 1);
390  break;
391  case math::TVD_RK3:
392  // Perform one explicit Euler step: t1 = t0 + dt
393  // Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]|
394  mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, /*result=*/1, /*speed*/3);
395  // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
396  this->cook(PARALLEL_FOR, 1);
397 
398  // Convex combine explict Euler step: t2 = t0 + dt/2
399  // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|)
400  mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.75),
401  /*phi=*/1, /*result=*/2, /*speed*/3);
402  // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
403  this->cook(PARALLEL_FOR, 2);
404 
405  // Convex combine explict Euler step: t3 = t0 + dt
406  // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|)
407  mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(1.0/3.0),
408  /*phi=*/1, /*result=*/2, /*speed*/3);
409  // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
410  this->cook(PARALLEL_FOR, 2);
411  break;
412  default:
413  OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
414  }//end of compile-time resolved switch
416 
417  time0 += dt;
418  ++countCFL;
419  mTracker->leafs().removeAuxBuffers();
420 
421  // Track the narrow band
422  mTracker->track();
423  }//end wile-loop over time
424 
425  return countCFL;//number of CLF propagation steps
426 }
427 
428 template<typename GridT, typename InterruptT>
429 template<typename MapT, math::BiasedGradientScheme SpatialScheme,
430  math::TemporalIntegrationScheme TemporalScheme>
431 inline typename GridT::ValueType
432 LevelSetMorphing<GridT, InterruptT>::
433 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
434 sampleSpeed(ScalarType time0, ScalarType time1, Index speedBuffer)
435 {
436  mMaxAbsS = mMinAbsS;
437  const size_t leafCount = mTracker->leafs().leafCount();
438  if (leafCount==0 || time0 >= time1) return ScalarType(0);
439 
440  if (mTarget->transform() == mTracker->grid().transform()) {
441  mTask = boost::bind(&LevelSetMorph::sampleAlignedSpeed, _1, _2, speedBuffer);
442  } else {
443  mTask = boost::bind(&LevelSetMorph::sampleXformedSpeed, _1, _2, speedBuffer);
444  }
445  this->cook(PARALLEL_REDUCE);
446  if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ScalarType(0);//speed is essentially zero
447  static const ScalarType CFL = (TemporalScheme == math::TVD_RK1 ? ScalarType(0.3) :
448  TemporalScheme == math::TVD_RK2 ? ScalarType(0.9) :
449  ScalarType(1.0))/math::Sqrt(ScalarType(3.0));
450  const ScalarType dt = math::Abs(time1 - time0), dx = mTracker->voxelSize();
451  return math::Min(dt, ScalarType(CFL*dx/mMaxAbsS));
452 }
453 
454 template<typename GridT, typename InterruptT>
455 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
456  math::TemporalIntegrationScheme TemporalScheme>
457 inline void
458 LevelSetMorphing<GridT, InterruptT>::
459 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
460 sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
461 {
462  typedef typename LeafType::ValueOnCIter VoxelIterT;
463  typedef tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler> SamplerT;
464  const MapT& map = *mMap;
465  mTracker->checkInterrupter();
466 
467  SamplerT sampler(mTarget->getAccessor(), mTarget->transform());
468  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
469  BufferType& speed = leafIter.buffer(speedBuffer);
470  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
471  ScalarType& s = const_cast<ScalarType&>(speed.getValue(voxelIter.pos()));
472  s -= sampler.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
473  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
474  }
475  }
476 }
477 
478 template<typename GridT, typename InterruptT>
479 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
480  math::TemporalIntegrationScheme TemporalScheme>
481 inline void
482 LevelSetMorphing<GridT, InterruptT>::
483 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
484 sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
485 {
486  typedef typename LeafType::ValueOnCIter VoxelIterT;
487  mTracker->checkInterrupter();
488 
489  typename GridT::ConstAccessor target = mTarget->getAccessor();
490  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
491  BufferType& speed = leafIter.buffer(speedBuffer);
492  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
493  ScalarType& s = const_cast<ScalarType&>(speed.getValue(voxelIter.pos()));
494  s -= target.getValue(voxelIter.getCoord());
495  mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
496  }
497  }
498 }
499 
500 template<typename GridT, typename InterruptT>
501 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
502  math::TemporalIntegrationScheme TemporalScheme>
503 inline void
504 LevelSetMorphing<GridT, InterruptT>::
505 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
506 cook(ThreadingMode mode, size_t swapBuffer)
507 {
508  mTracker->startInterrupter("Morphing level set");
509 
510  const int grainSize = mTracker->getGrainSize();
511  const LeafRange range = mTracker->leafs().leafRange(grainSize);
512 
513  if (mTracker->getGrainSize()==0) {
514  (*this)(range);
515  } else if (mode == PARALLEL_FOR) {
516  tbb::parallel_for(range, *this);
517  } else if (mode == PARALLEL_REDUCE) {
518  tbb::parallel_reduce(range, *this);
519  } else {
520  throw std::runtime_error("Undefined threading mode");
521  }
522 
523  mTracker->leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
524 
525  mTracker->endInterrupter();
526 }
527 
528 // Forward Euler advection steps:
529 // Phi(result) = Phi(0) - dt * Phi(speed) * |Grad[Phi(0)]|
530 template<typename GridT, typename InterruptT>
531 template <typename MapT,
532  math::BiasedGradientScheme SpatialScheme,
533  math::TemporalIntegrationScheme TemporalScheme>
534 inline void
535 LevelSetMorphing<GridT, InterruptT>::
536 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
537 euler1(const LeafRange& range, ScalarType dt, Index resultBuffer, Index speedBuffer)
538 {
539  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
540  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
541  typedef typename LeafType::ValueOnCIter VoxelIterT;
542 
543  mTracker->checkInterrupter();
544  const MapT& map = *mMap;
545  Stencil stencil(mTracker->grid());
546 
547  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
548  BufferType& speed = leafIter.buffer(speedBuffer);
549  BufferType& result = leafIter.buffer(resultBuffer);
550  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
551  const Index n = voxelIter.pos();
552  stencil.moveTo(voxelIter);
553  const ScalarType G = math::GradientNormSqrd<MapT,SpatialScheme>::result(map, stencil);
554  result.setValue(n, *voxelIter - dt * speed.getValue(n) * G);
555  }
556  }
557 }
558 
559 // Convex combination of Phi and a forward Euler advection steps:
560 // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Phi(speed) * |Grad[Phi(0)]|)
561 template<typename GridT, typename InterruptT>
562 template <typename MapT, math::BiasedGradientScheme SpatialScheme,
563  math::TemporalIntegrationScheme TemporalScheme>
564 inline void
565 LevelSetMorphing<GridT, InterruptT>::
566 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
567 euler2(const LeafRange& range, ScalarType dt, ScalarType alpha,
568  Index phiBuffer, Index resultBuffer, Index speedBuffer)
569 {
570  typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
571  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
572  typedef typename LeafType::ValueOnCIter VoxelIterT;
573 
574  mTracker->checkInterrupter();
575  const MapT& map = *mMap;
576  const ScalarType beta = ScalarType(1.0) - alpha;
577  Stencil stencil(mTracker->grid());
578 
579  for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
580  BufferType& speed = leafIter.buffer(speedBuffer);
581  BufferType& result = leafIter.buffer(resultBuffer);
582  BufferType& phi = leafIter.buffer(phiBuffer);
583  for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
584  const Index n = voxelIter.pos();
585  stencil.moveTo(voxelIter);
586  const ScalarType G = math::GradientNormSqrd<MapT,SpatialScheme>::result(map, stencil);
587  result.setValue(n,
588  alpha * phi.getValue(n) + beta * (*voxelIter - dt * speed.getValue(n) * G));
589  }
590  }
591 }
592 
593 } // namespace tools
594 } // namespace OPENVDB_VERSION_NAME
595 } // namespace openvdb
596 
597 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
598 
599 // Copyright (c) 2012-2013 DreamWorks Animation LLC
600 // All rights reserved. This software is distributed under the
601 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
virtual ~LevelSetMorphing()
Definition: LevelSetMorph.h:91
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
void setGrainSize(int grainsize)
Set the grain size used for multithreading.
Definition: LevelSetMorph.h:135
TrackerT::BufferType BufferType
Definition: LevelSetMorph.h:82
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:848
TrackerT::LeafRange LeafRange
Definition: LevelSetMorph.h:80
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:653
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=NULL)
Main constructor.
Definition: LevelSetMorph.h:86
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:112
Index32 Index
Definition: Types.h:57
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Name mapType() const
Return the transformation map's type-name.
Definition: Transform.h:92
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:122
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:244
TrackerT::LeafType LeafType
Definition: LevelSetMorph.h:81
TrackerT::ValueType ScalarType
Definition: LevelSetMorph.h:83
Definition: FiniteDifference.h:264
LeafManagerType::BufferType BufferType
Definition: LevelSetTracker.h:76
Definition: Exceptions.h:88
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition: LevelSetMorph.h:104
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:229
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:102
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:505
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: FiniteDifference.h:265
GridT::TreeType TreeType
Definition: LevelSetMorph.h:78
#define OPENVDB_VERSION_NAME
Definition: version.h:45
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition: LevelSetMorph.h:99
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:129
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
TreeType::LeafNodeType LeafType
Definition: LevelSetTracker.h:71
Definition: FiniteDifference.h:198
void setTarget(const GridT &targetGrid)
Redefine the target level set.
Definition: LevelSetMorph.h:94
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:121
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Return the temporal integration scheme.
Definition: LevelSetMorph.h:117
GridT GridType
Definition: LevelSetMorph.h:77
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:107
TreeType::ValueType ValueType
Definition: LevelSetTracker.h:72
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
Definition: LevelSetMorph.h:97
int getNormCount() const
Return the number of normalizations performed per track or normalize call.
Definition: LevelSetMorph.h:127
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:122
LevelSetTracker< GridT, InterruptT > TrackerT
Definition: LevelSetMorph.h:79
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition: LevelSetMorph.h:74
Definition: FiniteDifference.h:263
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:66
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:274
int getGrainSize() const
Return the grain size used for multithreading.
Definition: LevelSetMorph.h:132
Definition: FiniteDifference.h:194
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:65
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:566