12 #ifndef OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
13 #define OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
15 #include <tbb/parallel_for.h>
30 #include <type_traits>
54 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
62 using LeafType =
typename TreeType::LeafNodeType;
65 using LeafRange =
typename LeafManagerType::LeafRange;
67 using MaskTreeType =
typename TreeType::template ValueConverter<ValueMask>::Type;
68 static_assert(std::is_floating_point<ValueType>::value,
69 "LevelSetTracker requires a level set grid with floating-point values");
76 : spatialScheme(s), temporalScheme(t), normCount(n), grainSize(g) {}
91 template <
typename MaskType>
95 void normalize() { this->normalize<MaskTreeType>(
nullptr); }
198 template<TrimMode Trimming>
203 void operator()(
const LeafRange& r)
const;
204 LevelSetTracker& mTracker;
214 using StencilT =
typename SchemeT::template ISStencil<GridType>::StencilType;
215 using MaskLeafT =
typename MaskT::LeafNodeType;
216 using MaskIterT =
typename MaskLeafT::ValueOnCIter;
217 using VoxelIterT =
typename LeafType::ValueOnCIter;
219 Normalizer(LevelSetTracker& tracker,
const MaskT* mask);
221 void operator()(
const LeafRange& r)
const {mTask(
const_cast<Normalizer*
>(
this), r);}
222 void cook(
const char* msg,
int swapBuffer=0);
223 template <
int Nominator,
int Denominator>
224 void euler(
const LeafRange& range,
Index phiBuffer,
Index resultBuffer);
225 inline void euler01(
const LeafRange& r) {this->euler<0,1>(r, 0, 1);}
226 inline void euler12(
const LeafRange& r) {this->euler<1,2>(r, 1, 1);}
227 inline void euler34(
const LeafRange& r) {this->euler<3,4>(r, 1, 2);}
228 inline void euler13(
const LeafRange& r) {this->euler<1,3>(r, 1, 2);}
229 template <
int Nominator,
int Denominator>
230 void eval(StencilT& stencil,
const ValueType* phi, ValueType* result,
Index n)
const;
231 LevelSetTracker& mTracker;
233 const ValueType mDt, mInvDx;
234 typename std::function<void (Normalizer*,
const LeafRange&)> mTask;
237 template<math::BiasedGradientScheme SpatialScheme,
typename MaskT>
238 void normalize1(
const MaskT* mask);
242 void normalize2(
const MaskT* mask);
249 LeafManagerType* mLeafs;
250 InterruptT* mInterrupter;
253 TrimMode mTrimMode = TrimMode::kAll;
256 template<
typename Gr
idT,
typename InterruptT>
261 mInterrupter(interrupt),
262 mDx(static_cast<
ValueType>(grid.voxelSize()[0])),
265 if ( !
grid.hasUniformVoxels() ) {
267 "The transform must have uniform scale for the LevelSetTracker to function");
271 "LevelSetTracker expected a level set, got a grid of class \""
272 +
grid.gridClassToString(
grid.getGridClass())
273 +
"\" [hint: Grid::setGridClass(openvdb::GRID_LEVEL_SET)]");
277 template<
typename Gr
idT,
typename InterruptT>
282 this->startInterrupter(
"Pruning Level Set");
286 case TrimMode::kNone:
break;
287 case TrimMode::kInterior: Trim<TrimMode::kInterior>(*this).trim();
break;
288 case TrimMode::kExterior: Trim<TrimMode::kExterior>(*this).trim();
break;
289 case TrimMode::kAll: Trim<TrimMode::kAll>(*this).trim();
break;
296 mLeafs->rebuildLeafArray();
297 this->endInterrupter();
300 template<
typename Gr
idT,
typename InterruptT>
315 template<
typename Gr
idT,
typename InterruptT>
320 if (this->getNormCount() == 0) {
321 for (
int i=0; i < iterations; ++i) {
326 for (
int i=0; i < iterations; ++i) {
331 mask.topologyDifference(mask0);
337 template<
typename Gr
idT,
typename InterruptT>
340 erode(
int iterations)
343 mLeafs->rebuildLeafArray();
348 template<
typename Gr
idT,
typename InterruptT>
353 const int wOld =
static_cast<int>(
math::RoundDown(this->getHalfWidth()));
354 const int wNew =
static_cast<int>(halfWidth);
356 this->dilate(wNew - wOld);
357 }
else if (wOld > wNew) {
358 this->erode(wOld - wNew);
363 template<
typename Gr
idT,
typename InterruptT>
368 if (mInterrupter) mInterrupter->start(msg);
371 template<
typename Gr
idT,
typename InterruptT>
376 if (mInterrupter) mInterrupter->end();
379 template<
typename Gr
idT,
typename InterruptT>
385 tbb::task::self().cancel_group_execution();
391 template<
typename Gr
idT,
typename InterruptT>
392 template<
typename MaskT>
397 switch (this->getSpatialScheme()) {
399 this->normalize1<math::FIRST_BIAS , MaskT>(mask);
break;
401 this->normalize1<math::SECOND_BIAS, MaskT>(mask);
break;
403 this->normalize1<math::THIRD_BIAS, MaskT>(mask);
break;
405 this->normalize1<math::WENO5_BIAS, MaskT>(mask);
break;
407 this->normalize1<math::HJWENO5_BIAS, MaskT>(mask);
break;
414 template<
typename Gr
idT,
typename InterruptT>
415 template<math::BiasedGradientScheme SpatialScheme,
typename MaskT>
420 switch (this->getTemporalScheme()) {
422 this->normalize2<SpatialScheme, math::TVD_RK1, MaskT>(mask);
break;
424 this->normalize2<SpatialScheme, math::TVD_RK2, MaskT>(mask);
break;
426 this->normalize2<SpatialScheme, math::TVD_RK3, MaskT>(mask);
break;
433 template<
typename Gr
idT,
typename InterruptT>
438 LevelSetTracker<GridT, InterruptT>::
439 normalize2(
const MaskT* mask)
441 Normalizer<SpatialScheme, TemporalScheme, MaskT> tmp(*
this, mask);
449 template<
typename Gr
idT,
typename InterruptT>
450 template<lstrack::TrimMode Trimming>
452 LevelSetTracker<GridT, InterruptT>::Trim<Trimming>::trim()
455 if (Trimming != TrimMode::kNone) {
456 const int grainSize = mTracker.getGrainSize();
457 const LeafRange range = mTracker.leafs().leafRange(grainSize);
460 tbb::parallel_for(range, *
this);
470 template<
typename Gr
idT,
typename InterruptT>
471 template<lstrack::TrimMode Trimming>
473 LevelSetTracker<GridT, InterruptT>::Trim<Trimming>::operator()(
const LeafRange& range)
const
475 mTracker.checkInterrupter();
476 const ValueType gamma = mTracker.mGrid->background();
479 for (
auto leafIter = range.begin(); leafIter; ++leafIter) {
480 auto& leaf = *leafIter;
481 for (
auto iter = leaf.beginValueOn(); iter; ++iter) {
482 const auto val = *iter;
484 case TrimMode::kNone:
486 case TrimMode::kInterior:
487 if (val <= -gamma) { leaf.setValueOff(iter.pos(), -gamma); }
489 case TrimMode::kExterior:
490 if (val >= gamma) { leaf.setValueOff(iter.pos(), gamma); }
494 leaf.setValueOff(iter.pos(), -gamma);
495 }
else if (val >= gamma) {
496 leaf.setValueOff(iter.pos(), gamma);
508 template<
typename Gr
idT,
typename InterruptT>
513 LevelSetTracker<GridT, InterruptT>::
514 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
515 Normalizer(LevelSetTracker& tracker,
const MaskT* mask)
518 , mDt(tracker.voxelSize()*(TemporalScheme == math::
TVD_RK1 ? 0.3f :
519 TemporalScheme == math::
TVD_RK2 ? 0.9f : 1.0f))
520 , mInvDx(1.0f/tracker.voxelSize())
525 template<
typename Gr
idT,
typename InterruptT>
534 namespace ph = std::placeholders;
537 mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
539 for (
int n=0, e=mTracker.getNormCount(); n < e; ++n) {
542 switch(TemporalScheme) {
546 mTask = std::bind(&Normalizer::euler01, ph::_1, ph::_2);
549 this->cook(
"Normalizing level set using TVD_RK1", 1);
554 mTask = std::bind(&Normalizer::euler01, ph::_1, ph::_2);
557 this->cook(
"Normalizing level set using TVD_RK1 (step 1 of 2)", 1);
561 mTask = std::bind(&Normalizer::euler12, ph::_1, ph::_2);
564 this->cook(
"Normalizing level set using TVD_RK1 (step 2 of 2)", 1);
569 mTask = std::bind(&Normalizer::euler01, ph::_1, ph::_2);
572 this->cook(
"Normalizing level set using TVD_RK3 (step 1 of 3)", 1);
576 mTask = std::bind(&Normalizer::euler34, ph::_1, ph::_2);
579 this->cook(
"Normalizing level set using TVD_RK3 (step 2 of 3)", 2);
583 mTask = std::bind(&Normalizer::euler13, ph::_1, ph::_2);
586 this->cook(
"Normalizing level set using TVD_RK3 (step 3 of 3)", 2);
590 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
594 mTracker.mLeafs->removeAuxBuffers();
599 template<
typename Gr
idT,
typename InterruptT>
604 LevelSetTracker<GridT, InterruptT>::
605 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
606 cook(
const char* msg,
int swapBuffer)
608 mTracker.startInterrupter( msg );
610 const int grainSize = mTracker.getGrainSize();
611 const LeafRange range = mTracker.leafs().leafRange(grainSize);
613 grainSize>0 ? tbb::parallel_for(range, *
this) : (*this)(range);
615 mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize==0);
617 mTracker.endInterrupter();
620 template<
typename Gr
idT,
typename InterruptT>
624 template <
int Nominator,
int Denominator>
626 LevelSetTracker<GridT, InterruptT>::
627 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
628 eval(StencilT& stencil,
const ValueType* phi, ValueType* result,
Index n)
const
630 using GradientT =
typename math::ISGradientNormSqrd<SpatialScheme>;
631 static const ValueType alpha = ValueType(Nominator)/ValueType(Denominator);
632 static const ValueType beta = ValueType(1) - alpha;
634 const ValueType normSqGradPhi = GradientT::result(stencil);
635 const ValueType phi0 = stencil.getValue();
637 math::Tolerance<ValueType>::value() );
638 v = phi0 - mDt * v * (
math::Sqrt(normSqGradPhi) * mInvDx - 1.0f);
639 result[n] = Nominator ? alpha * phi[n] + beta * v : v;
642 template<
typename Gr
idT,
typename InterruptT>
646 template <
int Nominator,
int Denominator>
648 LevelSetTracker<GridT,InterruptT>::
649 Normalizer<SpatialScheme, TemporalScheme, MaskT>::
650 euler(
const LeafRange& range,
Index phiBuffer,
Index resultBuffer)
652 mTracker.checkInterrupter();
654 StencilT stencil(mTracker.grid());
656 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
657 const ValueType* phi = leafIter.buffer(phiBuffer).data();
658 ValueType* result = leafIter.buffer(resultBuffer).data();
659 if (mMask ==
nullptr) {
660 for (
auto iter = leafIter->cbeginValueOn(); iter; ++iter) {
661 stencil.moveTo(iter);
662 this->eval<Nominator, Denominator>(stencil, phi, result, iter.pos());
664 }
else if (
const MaskLeafT* mask = mMask->probeLeaf(leafIter->origin())) {
665 const ValueType* phi0 = leafIter->buffer().data();
666 for (MaskIterT iter = mask->cbeginValueOn(); iter; ++iter) {
667 const Index i = iter.pos();
668 stencil.moveTo(iter.getCoord(), phi0[i]);
669 this->eval<Nominator, Denominator>(stencil, phi, result, i);
679 #endif // OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED