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) {
370 switch(TemporalScheme) {
375 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
383 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
389 mTask = boost::bind(&Normalizer::euler2,
390 _1, _2, dt, ValueType(0.5), 1, 1);
398 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
404 mTask = boost::bind(&Normalizer::euler2,
405 _1, _2, dt, ValueType(0.75), 1, 2);
411 mTask = boost::bind(&Normalizer::euler2,
412 _1, _2, dt, ValueType(1.0/3.0), 1, 2);
417 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
421 mTracker.mLeafs->removeAuxBuffers();
426 template<
typename Gr
idT,
typename InterruptT>
430 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
433 mTracker.startInterrupter(
"Normalizing Level Set");
435 if (mTracker.getGrainSize()>0) {
436 tbb::parallel_for(mTracker.mLeafs->getRange(mTracker.getGrainSize()), *
this);
438 (*this)(mTracker.mLeafs->getRange());
441 mTracker.mLeafs->swapLeafBuffer(swapBuffer, mTracker.getGrainSize()==0);
443 mTracker.endInterrupter();
450 template<
typename Gr
idT,
typename InterruptT>
454 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
455 euler1(
const RangeType &range, ValueType dt,
Index resultBuffer)
457 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
458 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
459 typedef typename LeafType::ValueOnCIter VoxelIterT;
460 mTracker.checkInterrupter();
461 const ValueType one(1.0), invDx = one/mTracker.voxelSize();
462 Stencil stencil(mTracker.grid());
463 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
464 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
465 const LeafType& leaf = mTracker.mLeafs->leaf(n);
466 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
467 stencil.moveTo(iter);
468 const ValueType normSqGradPhi =
470 const ValueType phi0 = stencil.getValue();
471 const ValueType diff =
math::Sqrt(normSqGradPhi)*invDx - one;
473 result.setValue(iter.pos(), phi0 - dt * S * diff);
478 template<
typename Gr
idT,
typename InterruptT>
482 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
483 euler2(
const RangeType& range, ValueType dt, ValueType alpha,
Index phiBuffer,
Index resultBuffer)
485 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
486 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
487 typedef typename LeafType::ValueOnCIter VoxelIterT;
488 mTracker.checkInterrupter();
489 const ValueType one(1.0), beta = one - alpha, invDx = one/mTracker.voxelSize();
490 Stencil stencil(mTracker.grid());
491 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
492 const BufferType& phi = mTracker.mLeafs->getBuffer(n, phiBuffer);
493 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
494 const LeafType& leaf = mTracker.mLeafs->leaf(n);
495 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
496 stencil.moveTo(iter);
497 const ValueType normSqGradPhi =
499 const ValueType phi0 = stencil.getValue();
500 const ValueType diff =
math::Sqrt(normSqGradPhi)*invDx - one;
502 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(phi0 - dt * S * diff));
511 #endif // OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
Definition: FiniteDifference.h:264
Definition: FiniteDifference.h:263
Definition: FiniteDifference.h:195
TemporalIntegrationScheme
Temporal integrations schemes.
Definition: FiniteDifference.h:261
Definition: FiniteDifference.h:198
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:648
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Definition: FiniteDifference.h:194
Type Pow2(Type x)
Return .
Definition: Math.h:460
Definition: Exceptions.h:86
Definition: FiniteDifference.h:265
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:260
#define OPENVDB_VERSION_NAME
Definition: version.h:45
tbb::blocked_range< size_t > RangeType
Definition: LeafManager.h:118
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:117
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:143
Index32 Index
Definition: Types.h:56
Definition: FiniteDifference.h:197
Definition: Exceptions.h:88
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Definition: FiniteDifference.h:196
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:56
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76