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>
55 #include <openvdb/tree/ValueAccessor.h>
56 #include <openvdb/tree/LeafManager.h>
65 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
71 typedef typename TreeType::LeafNodeType
LeafType;
77 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
123 void startInterrupter(
const char* msg);
124 void endInterrupter();
126 bool checkInterrupter();
137 if (mTask) mTask(const_cast<LevelSetTracker*>(
this), r);
149 void operator()(
const RangeType& r)
const {mTask(const_cast<Normalizer*>(
this), r);}
150 typedef typename boost::function<void (Normalizer*, const RangeType&)> FuncType;
151 LevelSetTracker& mTracker;
153 void cook(
int swapBuffer=0);
154 void euler1(
const RangeType& range, ValueType dt,
Index resultBuffer);
155 void euler2(
const RangeType& range, ValueType dt, ValueType alpha,
159 typedef typename boost::function<void (LevelSetTracker*, const RangeType&)> FuncType;
161 void trim(
const RangeType& r);
163 template<math::BiasedGradientScheme SpatialScheme>
175 LeafManagerType* mLeafs;
176 InterruptT* mInterrupter;
183 const bool mIsMaster;
186 void operator=(
const LevelSetTracker& other) {}
190 template<
typename Gr
idT,
typename InterruptT>
194 mInterrupter(interrupt),
195 mDx(grid.voxelSize()[0]),
197 mTemporalScheme(math::
TVD_RK1),
203 if ( !grid.hasUniformVoxels() ) {
205 "The transform must have uniform scale for the LevelSetTracker to function");
209 "LevelSetTracker only supports level sets!\n"
210 "However, only level sets are guaranteed to work!\n"
211 "Hint: Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
215 template<
typename Gr
idT,
typename InterruptT>
218 mLeafs(other.mLeafs),
219 mInterrupter(other.mInterrupter),
221 mSpatialScheme(other.mSpatialScheme),
222 mTemporalScheme(other.mTemporalScheme),
223 mNormCount(other.mNormCount),
224 mGrainSize(other.mGrainSize),
230 template<
typename Gr
idT,
typename InterruptT>
234 this->startInterrupter(
"Pruning Level Set");
236 mTask = boost::bind(&LevelSetTracker::trim, _1, _2);
238 tbb::parallel_for(mLeafs->getRange(mGrainSize), *
this);
240 (*this)(mLeafs->getRange());
244 mGrid->tree().pruneLevelSet();
247 mLeafs->rebuildLeafArray();
248 this->endInterrupter();
251 template<
typename Gr
idT,
typename InterruptT>
265 template<
typename Gr
idT,
typename InterruptT>
269 if (mInterrupter) mInterrupter->start(msg);
272 template<
typename Gr
idT,
typename InterruptT>
276 if (mInterrupter) mInterrupter->end();
279 template<
typename Gr
idT,
typename InterruptT>
284 tbb::task::self().cancel_group_execution();
291 template<
typename Gr
idT,
typename InterruptT>
295 typedef typename LeafType::ValueOnIter VoxelIterT;
297 const ValueType gamma = mGrid->background();
298 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
299 LeafType &leaf = mLeafs->leaf(n);
300 for (VoxelIterT iter = leaf.beginValueOn(); iter; ++iter) {
301 const ValueType val = *iter;
303 leaf.setValueOff(iter.pos(), -gamma);
304 else if (val > gamma)
305 leaf.setValueOff(iter.pos(), gamma);
310 template<
typename Gr
idT,
typename InterruptT>
314 switch (mSpatialScheme) {
316 this->normalize1<math::FIRST_BIAS >();
break;
318 this->normalize1<math::SECOND_BIAS >();
break;
320 this->normalize1<math::THIRD_BIAS >();
break;
322 this->normalize1<math::WENO5_BIAS >();
break;
324 this->normalize1<math::HJWENO5_BIAS>();
break;
330 template<
typename Gr
idT,
typename InterruptT>
331 template<math::BiasedGradientScheme SpatialScheme>
335 switch (mTemporalScheme) {
337 this->normalize2<SpatialScheme, math::TVD_RK1>();
break;
339 this->normalize2<SpatialScheme, math::TVD_RK2>();
break;
341 this->normalize2<SpatialScheme, math::TVD_RK3>();
break;
347 template<
typename Gr
idT,
typename InterruptT>
351 LevelSetTracker<GridT, InterruptT>::normalize2()
353 Normalizer<SpatialScheme, TemporalScheme> tmp(*
this);
357 template<
typename Gr
idT,
typename InterruptT>
365 mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
367 const ValueType dt = (TemporalScheme ==
math::TVD_RK1 ? ValueType(0.3) :
368 TemporalScheme == math::
TVD_RK2 ? ValueType(0.9) : ValueType(1.0))
369 * ValueType(mTracker.voxelSize());
371 for (
int n=0, e=mTracker.getNormCount(); n < e; ++n) {
374 switch(TemporalScheme) {
379 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
387 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
393 mTask = boost::bind(&Normalizer::euler2,
394 _1, _2, dt, ValueType(0.5), 1, 1);
402 mTask = boost::bind(&Normalizer::euler1, _1, _2, dt, 1);
408 mTask = boost::bind(&Normalizer::euler2,
409 _1, _2, dt, ValueType(0.75), 1, 2);
415 mTask = boost::bind(&Normalizer::euler2,
416 _1, _2, dt, ValueType(1.0/3.0), 1, 2);
421 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
425 mTracker.mLeafs->removeAuxBuffers();
430 template<
typename Gr
idT,
typename InterruptT>
434 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
437 mTracker.startInterrupter(
"Normalizing Level Set");
439 if (mTracker.getGrainSize()>0) {
440 tbb::parallel_for(mTracker.mLeafs->getRange(mTracker.getGrainSize()), *
this);
442 (*this)(mTracker.mLeafs->getRange());
445 mTracker.mLeafs->swapLeafBuffer(swapBuffer, mTracker.getGrainSize()==0);
447 mTracker.endInterrupter();
454 template<
typename Gr
idT,
typename InterruptT>
458 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
459 euler1(
const RangeType &range, ValueType dt,
Index resultBuffer)
461 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
462 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
463 typedef typename LeafType::ValueOnCIter VoxelIterT;
464 mTracker.checkInterrupter();
465 const ValueType one(1.0), invDx = one/mTracker.voxelSize();
466 Stencil stencil(mTracker.grid());
467 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
468 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
469 const LeafType& leaf = mTracker.mLeafs->leaf(n);
470 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
471 stencil.moveTo(iter);
472 const ValueType normSqGradPhi =
474 const ValueType phi0 = stencil.getValue();
475 const ValueType diff =
math::Sqrt(normSqGradPhi)*invDx - one;
477 result.setValue(iter.pos(), phi0 - dt * S * diff);
482 template<
typename Gr
idT,
typename InterruptT>
486 LevelSetTracker<GridT,InterruptT>::Normalizer<SpatialScheme, TemporalScheme>::
487 euler2(
const RangeType& range, ValueType dt, ValueType alpha,
Index phiBuffer,
Index resultBuffer)
489 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
490 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
491 typedef typename LeafType::ValueOnCIter VoxelIterT;
492 mTracker.checkInterrupter();
493 const ValueType one(1.0), beta = one - alpha, invDx = one/mTracker.voxelSize();
494 Stencil stencil(mTracker.grid());
495 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
496 const BufferType& phi = mTracker.mLeafs->getBuffer(n, phiBuffer);
497 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
498 const LeafType& leaf = mTracker.mLeafs->leaf(n);
499 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
500 stencil.moveTo(iter);
501 const ValueType normSqGradPhi =
503 const ValueType phi0 = stencil.getValue();
504 const ValueType diff =
math::Sqrt(normSqGradPhi)*invDx - one;
506 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(phi0 - dt * S * diff));
515 #endif // OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:653
Index32 Index
Definition: Types.h:57
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Definition: FiniteDifference.h:264
tbb::blocked_range< size_t > RangeType
Definition: LeafManager.h:118
Type Pow2(Type x)
Return .
Definition: Math.h:458
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:176
Definition: Exceptions.h:88
Definition: FiniteDifference.h:196
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: LeafManager.h:122
Definition: FiniteDifference.h:265
#define OPENVDB_VERSION_NAME
Definition: version.h:45
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:260
Definition: FiniteDifference.h:198
Definition: FiniteDifference.h:195
Definition: Exceptions.h:86
Definition: FiniteDifference.h:263
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:117
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:194