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>
67 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
73 typedef typename TreeType::LeafNodeType
LeafType;
79 typedef typename TreeType::template ValueConverter<bool>::Type
BoolMaskType;
80 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
93 template <
typename MaskType>
97 void normalize() { this->normalize<BoolMaskType>(NULL); }
108 void dilate(
int iterations = 1);
112 void erode(
int iterations = 1);
143 void startInterrupter(
const char* msg);
145 void endInterrupter();
148 bool checkInterrupter();
150 const GridType&
grid()
const {
return *mGrid; }
152 LeafManagerType&
leafs() {
return *mLeafs; }
154 const LeafManagerType&
leafs()
const {
return *mLeafs; }
164 void operator()(
const RangeType& r)
const;
165 LevelSetTracker& mTracker;
175 typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
176 typedef typename MaskT::LeafNodeType MaskLeafT;
177 typedef typename MaskLeafT::ValueOnCIter MaskIterT;
178 typedef typename LeafType::ValueOnCIter VoxelIterT;
179 Normalizer(LevelSetTracker& tracker,
const MaskT* mask);
181 void operator()(
const RangeType& r)
const {mTask(const_cast<Normalizer*>(
this), r);}
182 void cook(
int swapBuffer=0);
183 template <
int Nominator,
int Denominator>
184 void euler(
const RangeType& range,
Index phiBuffer,
Index resultBuffer);
185 inline void euler01(
const RangeType& r) {this->euler<0,1>(r, 0, 1);}
186 inline void euler12(
const RangeType& r,
Index n,
Index m) {this->euler<1,2>(r, n, m);}
187 inline void euler34(
const RangeType& r,
Index n,
Index m) {this->euler<3,4>(r, n, m);}
188 inline void euler13(
const RangeType& r,
Index n,
Index m) {this->euler<1,3>(r, n, m);}
189 template <
int Nominator,
int Denominator>
190 void eval(StencilT& stencil,
const BufferType& phi, BufferType& result,
Index n)
const;
191 typedef typename boost::function<void (Normalizer*, const RangeType&)> FuncType;
192 LevelSetTracker& mTracker;
194 const ValueType mDt, mInvDx;
198 template<math::BiasedGradientScheme SpatialScheme,
typename MaskT>
199 void normalize1(
const MaskT* mask);
203 void normalize2(
const MaskT* mask);
210 LeafManagerType* mLeafs;
211 InterruptT* mInterrupter;
217 const bool mIsMaster;
220 void operator=(
const LevelSetTracker& other) {}
224 template<
typename Gr
idT,
typename InterruptT>
225 LevelSetTracker<GridT, InterruptT>::
226 LevelSetTracker(GridT& grid, InterruptT* interrupt):
228 mLeafs(new LeafManagerType(grid.tree())),
229 mInterrupter(interrupt),
230 mDx(static_cast<ValueType>(grid.voxelSize()[0])),
232 mTemporalScheme(math::
TVD_RK1),
237 if ( !grid.hasUniformVoxels() ) {
239 "The transform must have uniform scale for the LevelSetTracker to function");
243 "LevelSetTracker expected a level set, got a grid of class \""
244 + grid.gridClassToString(grid.getGridClass())
245 +
"\" [hint: Grid::setGridClass(openvdb::GRID_LEVEL_SET)]");
249 template<
typename Gr
idT,
typename InterruptT>
253 mLeafs(other.mLeafs),
254 mInterrupter(other.mInterrupter),
256 mSpatialScheme(other.mSpatialScheme),
257 mTemporalScheme(other.mTemporalScheme),
258 mNormCount(other.mNormCount),
259 mGrainSize(other.mGrainSize),
264 template<
typename Gr
idT,
typename InterruptT>
269 this->startInterrupter(
"Pruning Level Set");
279 mLeafs->rebuildLeafArray();
280 this->endInterrupter();
283 template<
typename Gr
idT,
typename InterruptT>
298 template<
typename Gr
idT,
typename InterruptT>
303 if (mNormCount == 0) {
304 for (
int i=0; i < iterations; ++i) {
306 mLeafs->rebuildLeafArray();
310 for (
int i=0; i < iterations; ++i) {
311 BoolMaskType mask0(mGrid->tree(),
false,
TopologyCopy());
313 mLeafs->rebuildLeafArray();
315 BoolMaskType mask(mGrid->tree(),
false,
TopologyCopy());
316 mask.topologyDifference(mask0);
322 template<
typename Gr
idT,
typename InterruptT>
328 mLeafs->rebuildLeafArray();
329 const ValueType background = mGrid->background() - iterations*mDx;
333 template<
typename Gr
idT,
typename InterruptT>
338 if (mInterrupter) mInterrupter->start(msg);
341 template<
typename Gr
idT,
typename InterruptT>
346 if (mInterrupter) mInterrupter->end();
349 template<
typename Gr
idT,
typename InterruptT>
355 tbb::task::self().cancel_group_execution();
361 template<
typename Gr
idT,
typename InterruptT>
362 template<
typename MaskT>
367 switch (mSpatialScheme) {
369 this->normalize1<math::FIRST_BIAS , MaskT>(mask);
break;
371 this->normalize1<math::SECOND_BIAS, MaskT>(mask);
break;
373 this->normalize1<math::THIRD_BIAS, MaskT>(mask);
break;
375 this->normalize1<math::WENO5_BIAS, MaskT>(mask);
break;
377 this->normalize1<math::HJWENO5_BIAS, MaskT>(mask);
break;
383 template<
typename Gr
idT,
typename InterruptT>
384 template<math::BiasedGradientScheme SpatialScheme,
typename MaskT>
389 switch (mTemporalScheme) {
391 this->normalize2<SpatialScheme, math::TVD_RK1, MaskT>(mask);
break;
393 this->normalize2<SpatialScheme, math::TVD_RK2, MaskT>(mask);
break;
395 this->normalize2<SpatialScheme, math::TVD_RK3, MaskT>(mask);
break;
401 template<
typename Gr
idT,
typename InterruptT>
406 LevelSetTracker<GridT, InterruptT>::
407 normalize2(
const MaskT* mask)
409 Normalizer<SpatialScheme, TemporalScheme, MaskT> tmp(*
this, mask);
415 template<
typename Gr
idT,
typename InterruptT>
417 LevelSetTracker<GridT, InterruptT>::
420 const int grainSize = mTracker.getGrainSize();
422 tbb::parallel_for(mTracker.mLeafs->getRange(grainSize), *
this);
424 (*this)(mTracker.mLeafs->getRange());
429 template<
typename Gr
idT,
typename InterruptT>
431 LevelSetTracker<GridT, InterruptT>::
432 Trim::operator()(
const RangeType& range)
const
434 typedef typename LeafType::ValueOnIter VoxelIterT;
435 mTracker.checkInterrupter();
436 const ValueType gamma = mTracker.mGrid->background();
437 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
438 LeafType &leaf = mTracker.mLeafs->leaf(n);
439 for (VoxelIterT iter = leaf.beginValueOn(); iter; ++iter) {
440 const ValueType val = *iter;
442 leaf.setValueOff(iter.pos(), -gamma);
443 else if (val >= gamma)
444 leaf.setValueOff(iter.pos(), gamma);
451 template<
typename Gr
idT,
typename InterruptT>
456 LevelSetTracker<GridT, InterruptT>::
457 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
458 Normalizer(LevelSetTracker& tracker,
const MaskT* mask)
461 , mDt(tracker.voxelSize()*(TemporalScheme == math::
TVD_RK1 ? 0.3f :
462 TemporalScheme == math::
TVD_RK2 ? 0.9f : 1.0f))
463 , mInvDx(1.0f/tracker.voxelSize())
468 template<
typename Gr
idT,
typename InterruptT>
478 mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
480 for (
int n=0, e=mTracker.getNormCount(); n < e; ++n) {
483 switch(TemporalScheme) {
487 mTask = boost::bind(&Normalizer::euler01, _1, _2);
495 mTask = boost::bind(&Normalizer::euler01, _1, _2);
502 mTask = boost::bind(&Normalizer::euler12, _1, _2, 1, 1);
510 mTask = boost::bind(&Normalizer::euler01, _1, _2);
517 mTask = boost::bind(&Normalizer::euler34, _1, _2, 1, 2);
524 mTask = boost::bind(&Normalizer::euler13, _1, _2, 1, 2);
530 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
534 mTracker.mLeafs->removeAuxBuffers();
539 template<
typename Gr
idT,
typename InterruptT>
544 LevelSetTracker<GridT, InterruptT>::
545 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
548 mTracker.startInterrupter(
"Normalizing Level Set");
550 if (mTracker.getGrainSize()>0) {
551 tbb::parallel_for(mTracker.mLeafs->getRange(mTracker.getGrainSize()), *
this);
553 (*this)(mTracker.mLeafs->getRange());
556 mTracker.mLeafs->swapLeafBuffer(swapBuffer, mTracker.getGrainSize()==0);
558 mTracker.endInterrupter();
561 template<
typename Gr
idT,
typename InterruptT>
565 template <
int Nominator,
int Denominator>
567 LevelSetTracker<GridT, InterruptT>::
568 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
569 eval(StencilT& stencil,
const BufferType& phi, BufferType& result,
Index n)
const
571 typedef typename math::ISGradientNormSqrd<SpatialScheme> GradientT;
572 static const ValueType alpha = ValueType(Nominator)/ValueType(Denominator);
573 static const ValueType beta = ValueType(1) - alpha;
575 const ValueType normSqGradPhi = GradientT::result(stencil);
576 const ValueType phi0 = stencil.getValue();
578 v = phi0 - mDt * v * (
math::Sqrt(normSqGradPhi) * mInvDx - 1.0f);
579 result.setValue(n, Nominator ? alpha * phi[n] + beta * v : v);
582 template<
typename Gr
idT,
typename InterruptT>
586 template <
int Nominator,
int Denominator>
588 LevelSetTracker<GridT,InterruptT>::
589 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
590 euler(
const RangeType& range,
Index phiBuffer,
Index resultBuffer)
592 typedef typename LeafType::ValueOnCIter VoxelIterT;
594 mTracker.checkInterrupter();
596 StencilT stencil(mTracker.grid());
598 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
599 const BufferType& phi = mTracker.mLeafs->getBuffer(n, phiBuffer);
600 BufferType& result = mTracker.mLeafs->getBuffer(n, resultBuffer);
601 const LeafType& leaf = mTracker.mLeafs->leaf(n);
604 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter) {
605 stencil.moveTo(iter);
606 this->eval<Nominator,Denominator>(stencil, phi, result, iter.pos());
608 }
else if (
const MaskLeafT* mask = mMask->probeLeaf(leaf.origin())) {
609 for (MaskIterT iter = mask->cbeginValueOn(); iter; ++iter) {
610 const Index i = iter.pos();
611 stencil.moveTo(iter.getCoord(), leaf.getValue(i));
612 this->eval<Nominator,Denominator>(stencil, phi, result, i);
622 #endif // OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
Definition: LeafManager.h:126
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:265
tbb::blocked_range< size_t > RangeType
Definition: LeafManager.h:121
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Efficient multi-threaded replacement of the background values in tree.
Type Pow2(Type x)
Return .
Definition: Math.h:498
Definition: Operators.h:152
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:196
Defined various multi-threaded utility functions for trees.
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Definition: Exceptions.h:86
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:693
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Index32 Index
Definition: Types.h:59
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:216
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:39
Definition: FiniteDifference.h:195
Definition: FiniteDifference.h:264
Definition: Exceptions.h:88
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Definition: FiniteDifference.h:198
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Definition: FiniteDifference.h:194
Definition: FiniteDifference.h:263
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:120
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261