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