39 #ifndef OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
40 #define OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
42 #include <tbb/parallel_reduce.h>
43 #include <tbb/parallel_for.h>
44 #include <boost/bind.hpp>
45 #include <boost/function.hpp>
46 #include <boost/type_traits/is_floating_point.hpp>
47 #include <openvdb/Types.h>
48 #include <openvdb/math/Math.h>
49 #include <openvdb/math/FiniteDifference.h>
50 #include <openvdb/math/Operators.h>
51 #include <openvdb/math/Stencils.h>
52 #include <openvdb/math/Transform.h>
53 #include <openvdb/Grid.h>
54 #include <openvdb/util/NullInterrupter.h>
63 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
69 typedef typename TreeType::LeafNodeType
LeafType;
74 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
119 void startInterrupter(
const char* msg);
120 void endInterrupter();
122 bool checkInterrupter();
133 if (mTask) mTask(const_cast<LevelSetTracker*>(
this), r);
145 void operator()(
const RangeType& r)
const {mTask(const_cast<Normalizer*>(
this), r);}
146 typedef typename boost::function<void (Normalizer*, const RangeType&)> FuncType;
147 LevelSetTracker& mTracker;
149 void cook(
int swapBuffer=0);
150 void euler1(
const RangeType& range, ValueType dt,
Index resultBuffer);
151 void euler2(
const RangeType& range, ValueType dt, ValueType alpha,
155 typedef typename boost::function<void (LevelSetTracker*, const RangeType&)> FuncType;
157 void trim(
const RangeType& r);
159 template<math::BiasedGradientScheme SpatialScheme>
171 LeafManagerType* mLeafs;
172 InterruptT* mInterrupter;
179 const bool mIsMaster;
182 void operator=(
const LevelSetTracker& other) {}
186 template<
typename Gr
idT,
typename InterruptT>
190 mInterrupter(interrupt),
191 mDx(grid.voxelSize()[0]),
193 mTemporalScheme(math::
TVD_RK1),
199 if ( !grid.hasUniformVoxels() ) {
201 "The transform must have uniform scale for the LevelSetTracker to function");
205 "LevelSetTracker only supports level sets!\n"
206 "However, only level sets are guaranteed to work!\n"
207 "Hint: Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
211 template<
typename Gr
idT,
typename InterruptT>
214 mLeafs(other.mLeafs),
215 mInterrupter(other.mInterrupter),
217 mSpatialScheme(other.mSpatialScheme),
218 mTemporalScheme(other.mTemporalScheme),
219 mNormCount(other.mNormCount),
220 mGrainSize(other.mGrainSize),
226 template<
typename Gr
idT,
typename InterruptT>
230 this->startInterrupter(
"Pruning Level Set");
232 mTask = boost::bind(&LevelSetTracker::trim, _1, _2);
234 tbb::parallel_for(mLeafs->getRange(mGrainSize), *
this);
236 (*this)(mLeafs->getRange());
240 mGrid->tree().pruneLevelSet();
243 mLeafs->rebuildLeafArray();
244 this->endInterrupter();
247 template<
typename Gr
idT,
typename InterruptT>
261 template<
typename Gr
idT,
typename InterruptT>
265 if (mInterrupter) mInterrupter->start(msg);
268 template<
typename Gr
idT,
typename InterruptT>
272 if (mInterrupter) mInterrupter->end();
275 template<
typename Gr
idT,
typename InterruptT>
280 tbb::task::self().cancel_group_execution();
287 template<
typename Gr
idT,
typename InterruptT>
291 typedef typename LeafType::ValueOnIter VoxelIterT;
293 const ValueType gamma = mGrid->background();
294 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
295 LeafType &leaf = mLeafs->leaf(n);
296 for (VoxelIterT iter = leaf.beginValueOn(); iter; ++iter) {
297 const ValueType val = *iter;
299 leaf.setValueOff(iter.pos(), -gamma);
300 else if (val > gamma)
301 leaf.setValueOff(iter.pos(), gamma);
306 template<
typename Gr
idT,
typename InterruptT>
310 switch (mSpatialScheme) {
312 this->normalize1<math::FIRST_BIAS >();
break;
314 this->normalize1<math::SECOND_BIAS >();
break;
316 this->normalize1<math::THIRD_BIAS >();
break;
318 this->normalize1<math::WENO5_BIAS >();
break;
320 this->normalize1<math::HJWENO5_BIAS>();
break;
326 template<
typename Gr
idT,
typename InterruptT>
327 template<math::BiasedGradientScheme SpatialScheme>
331 switch (mTemporalScheme) {
333 this->normalize2<SpatialScheme, math::TVD_RK1>();
break;
335 this->normalize2<SpatialScheme, math::TVD_RK2>();
break;
337 this->normalize2<SpatialScheme, math::TVD_RK3>();
break;
343 template<
typename Gr
idT,
typename InterruptT>
347 LevelSetTracker<GridT, InterruptT>::normalize2()
349 Normalizer<SpatialScheme, TemporalScheme> tmp(*
this);
353 template<
typename Gr
idT,
typename InterruptT>
361 mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
363 const ValueType dt = (TemporalScheme ==
math::TVD_RK1 ? ValueType(0.3) :
364 TemporalScheme == math::
TVD_RK2 ? ValueType(0.9) : ValueType(1.0))
365 * ValueType(mTracker.voxelSize());
367 for (
int n=0, e=mTracker.getNormCount(); n < e; ++n) {
369 switch(TemporalScheme) {
374 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
382 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
388 mTask = boost::bind(&Normalizer::euler2,
389 _1, _2, dt, ValueType(0.5), 1, 1);
397 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
403 mTask = boost::bind(&Normalizer::euler2,
404 _1, _2, dt, ValueType(0.75), 1, 2);
410 mTask = boost::bind(&Normalizer::euler2,
411 _1, _2, dt, ValueType(1.0/3.0), 1, 2);
416 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
419 mTracker.mLeafs->removeAuxBuffers();
424 template<
typename Gr
idT,
typename InterruptT>
428 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
431 mTracker.startInterrupter(
"Normalizing Level Set");
433 if (mTracker.getGrainSize()>0) {
434 tbb::parallel_for(mTracker.mLeafs->getRange(mTracker.getGrainSize()), *
this);
436 (*this)(mTracker.mLeafs->getRange());
439 mTracker.mLeafs->swapLeafBuffer(swapBuffer, mTracker.getGrainSize()==0);
441 mTracker.endInterrupter();
448 template<
typename Gr
idT,
typename InterruptT>
452 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
453 euler1(
const RangeType &range, ValueType dt,
Index resultBuffer)
455 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
456 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
457 typedef typename LeafType::ValueOnCIter VoxelIterT;
458 mTracker.checkInterrupter();
459 const ValueType one(1.0), invDx = one/mTracker.voxelSize();
460 Stencil stencil(mTracker.grid());
461 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
462 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
463 const LeafType& leaf = mTracker.mLeafs->leaf(n);
464 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
465 stencil.moveTo(iter);
466 const ValueType normSqGradPhi =
468 const ValueType phi0 = stencil.getValue();
469 const ValueType diff =
math::Sqrt(normSqGradPhi)*invDx - one;
471 result.setValue(iter.pos(), phi0 - dt * S * diff);
476 template<
typename Gr
idT,
typename InterruptT>
480 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
481 euler2(
const RangeType& range, ValueType dt, ValueType alpha,
Index phiBuffer,
Index resultBuffer)
483 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
484 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
485 typedef typename LeafType::ValueOnCIter VoxelIterT;
486 mTracker.checkInterrupter();
487 const ValueType one(1.0), beta = one - alpha, invDx = one/mTracker.voxelSize();
488 Stencil stencil(mTracker.grid());
489 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
490 const BufferType& phi = mTracker.mLeafs->getBuffer(n, phiBuffer);
491 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
492 const LeafType& leaf = mTracker.mLeafs->leaf(n);
493 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
494 stencil.moveTo(iter);
495 const ValueType normSqGradPhi =
497 const ValueType phi0 = stencil.getValue();
498 const ValueType diff =
math::Sqrt(normSqGradPhi)*invDx - one;
500 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(phi0 - dt * S * diff));
509 #endif // OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED