12#ifndef OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
13#define OPENVDB_TOOLS_LEVEL_SET_TRACKER_HAS_BEEN_INCLUDED
23#include <openvdb/thread/Threading.h>
30#include <tbb/parallel_for.h>
56template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
64 using LeafType =
typename TreeType::LeafNodeType;
67 using LeafRange =
typename LeafManagerType::LeafRange;
69 using MaskTreeType =
typename TreeType::template ValueConverter<ValueMask>::Type;
70 static_assert(std::is_floating_point<ValueType>::value,
71 "LevelSetTracker requires a level set grid with floating-point values");
77 int n =
static_cast<int>(LEVEL_SET_HALF_WIDTH),
int g = 1)
78 : spatialScheme(s), temporalScheme(t), normCount(n), grainSize(g) {}
93 template <
typename MaskType>
97 void normalize() { this->normalize<MaskTreeType>(
nullptr); }
200 template<TrimMode Trimming>
205 void operator()(
const LeafRange& r)
const;
206 LevelSetTracker& mTracker;
216 using StencilT =
typename SchemeT::template ISStencil<GridType>::StencilType;
217 using MaskLeafT =
typename MaskT::LeafNodeType;
218 using MaskIterT =
typename MaskLeafT::ValueOnCIter;
219 using VoxelIterT =
typename LeafType::ValueOnCIter;
221 Normalizer(LevelSetTracker& tracker,
const MaskT* mask);
223 void operator()(
const LeafRange& r)
const {mTask(
const_cast<Normalizer*
>(
this), r);}
224 void cook(
const char* msg,
int swapBuffer=0);
225 template <
int Nominator,
int Denominator>
226 void euler(
const LeafRange& range, Index phiBuffer, Index resultBuffer);
227 inline void euler01(
const LeafRange& r) {this->euler<0,1>(r, 0, 1);}
228 inline void euler12(
const LeafRange& r) {this->euler<1,2>(r, 1, 1);}
229 inline void euler34(
const LeafRange& r) {this->euler<3,4>(r, 1, 2);}
230 inline void euler13(
const LeafRange& r) {this->euler<1,3>(r, 1, 2);}
231 template <
int Nominator,
int Denominator>
232 void eval(StencilT& stencil,
const ValueType* phi, ValueType* result, Index n)
const;
233 LevelSetTracker& mTracker;
235 const ValueType mDt, mInvDx;
236 typename std::function<void (Normalizer*,
const LeafRange&)> mTask;
239 template<math::BiasedGradientScheme SpatialScheme,
typename MaskT>
240 void normalize1(
const MaskT* mask);
242 template<math::BiasedGradientScheme SpatialScheme,
243 math::TemporalIntegrationScheme TemporalScheme,
typename MaskT>
244 void normalize2(
const MaskT* mask);
251 LeafManagerType* mLeafs;
252 InterruptT* mInterrupter;
255 TrimMode mTrimMode = TrimMode::kAll;
258template<
typename Gr
idT,
typename InterruptT>
263 mInterrupter(interrupt),
264 mDx(static_cast<
ValueType>(grid.voxelSize()[0])),
267 if ( !
grid.hasUniformVoxels() ) {
268 OPENVDB_THROW(RuntimeError,
269 "The transform must have uniform scale for the LevelSetTracker to function");
272 OPENVDB_THROW(RuntimeError,
273 "LevelSetTracker expected a level set, got a grid of class \""
274 + grid.gridClassToString(grid.getGridClass())
275 +
"\" [hint: Grid::setGridClass(openvdb::GRID_LEVEL_SET)]");
279template<
typename Gr
idT,
typename InterruptT>
284 this->startInterrupter(
"Pruning Level Set");
288 case TrimMode::kNone:
break;
289 case TrimMode::kInterior: Trim<TrimMode::kInterior>(*this).trim();
break;
290 case TrimMode::kExterior: Trim<TrimMode::kExterior>(*this).trim();
break;
291 case TrimMode::kAll: Trim<TrimMode::kAll>(*this).trim();
break;
298 mLeafs->rebuildLeafArray();
299 this->endInterrupter();
302template<
typename Gr
idT,
typename InterruptT>
317template<
typename Gr
idT,
typename InterruptT>
322 if (this->getNormCount() == 0) {
323 for (
int i=0; i < iterations; ++i) {
328 for (
int i=0; i < iterations; ++i) {
333 mask.topologyDifference(mask0);
339template<
typename Gr
idT,
typename InterruptT>
346 mLeafs->rebuildLeafArray();
351template<
typename Gr
idT,
typename InterruptT>
356 const int wOld =
static_cast<int>(
math::RoundDown(this->getHalfWidth()));
357 const int wNew =
static_cast<int>(halfWidth);
359 this->dilate(wNew - wOld);
360 }
else if (wOld > wNew) {
361 this->erode(wOld - wNew);
366template<
typename Gr
idT,
typename InterruptT>
371 if (mInterrupter) mInterrupter->start(msg);
374template<
typename Gr
idT,
typename InterruptT>
379 if (mInterrupter) mInterrupter->end();
382template<
typename Gr
idT,
typename InterruptT>
388 thread::cancelGroupExecution();
394template<
typename Gr
idT,
typename InterruptT>
395template<
typename MaskT>
400 switch (this->getSpatialScheme()) {
402 this->normalize1<math::FIRST_BIAS , MaskT>(mask);
break;
404 this->normalize1<math::SECOND_BIAS, MaskT>(mask);
break;
406 this->normalize1<math::THIRD_BIAS, MaskT>(mask);
break;
408 this->normalize1<math::WENO5_BIAS, MaskT>(mask);
break;
410 this->normalize1<math::HJWENO5_BIAS, MaskT>(mask);
break;
417template<
typename Gr
idT,
typename InterruptT>
418template<math::BiasedGradientScheme SpatialScheme,
typename MaskT>
423 switch (this->getTemporalScheme()) {
425 this->normalize2<SpatialScheme, math::TVD_RK1, MaskT>(mask);
break;
427 this->normalize2<SpatialScheme, math::TVD_RK2, MaskT>(mask);
break;
429 this->normalize2<SpatialScheme, math::TVD_RK3, MaskT>(mask);
break;
436template<
typename Gr
idT,
typename InterruptT>
441LevelSetTracker<GridT, InterruptT>::
442normalize2(
const MaskT* mask)
444 Normalizer<SpatialScheme, TemporalScheme, MaskT> tmp(*
this, mask);
452template<
typename Gr
idT,
typename InterruptT>
453template<lstrack::TrimMode Trimming>
455LevelSetTracker<GridT, InterruptT>::Trim<Trimming>::trim()
458 if (Trimming != TrimMode::kNone) {
459 const int grainSize = mTracker.getGrainSize();
460 const LeafRange range = mTracker.leafs().leafRange(grainSize);
463 tbb::parallel_for(range, *
this);
473template<
typename Gr
idT,
typename InterruptT>
474template<lstrack::TrimMode Trimming>
476LevelSetTracker<GridT, InterruptT>::Trim<Trimming>::operator()(
const LeafRange& range)
const
478 mTracker.checkInterrupter();
479 const ValueType gamma = mTracker.mGrid->background();
482 for (
auto leafIter = range.begin(); leafIter; ++leafIter) {
483 auto& leaf = *leafIter;
484 for (
auto iter = leaf.beginValueOn(); iter; ++iter) {
485 const auto val = *iter;
487 case TrimMode::kNone:
489 case TrimMode::kInterior:
490 if (val <= -gamma) { leaf.setValueOff(iter.pos(), -gamma); }
492 case TrimMode::kExterior:
493 if (val >= gamma) { leaf.setValueOff(iter.pos(), gamma); }
497 leaf.setValueOff(iter.pos(), -gamma);
498 }
else if (val >= gamma) {
499 leaf.setValueOff(iter.pos(), gamma);
511template<
typename Gr
idT,
typename InterruptT>
516LevelSetTracker<GridT, InterruptT>::
517Normalizer<SpatialScheme, TemporalScheme, MaskT>::
518Normalizer(LevelSetTracker& tracker,
const MaskT* mask)
521 , mDt(tracker.voxelSize()*(TemporalScheme == math::
TVD_RK1 ? 0.3f :
522 TemporalScheme == math::
TVD_RK2 ? 0.9f : 1.0f))
523 , mInvDx(1.0f/tracker.voxelSize())
528template<
typename Gr
idT,
typename InterruptT>
537 namespace ph = std::placeholders;
540 mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
542 for (
int n=0, e=mTracker.getNormCount(); n < e; ++n) {
545 switch(TemporalScheme) {
549 mTask = std::bind(&Normalizer::euler01, ph::_1, ph::_2);
552 this->cook(
"Normalizing level set using TVD_RK1", 1);
557 mTask = std::bind(&Normalizer::euler01, ph::_1, ph::_2);
560 this->cook(
"Normalizing level set using TVD_RK1 (step 1 of 2)", 1);
564 mTask = std::bind(&Normalizer::euler12, ph::_1, ph::_2);
567 this->cook(
"Normalizing level set using TVD_RK1 (step 2 of 2)", 1);
572 mTask = std::bind(&Normalizer::euler01, ph::_1, ph::_2);
575 this->cook(
"Normalizing level set using TVD_RK3 (step 1 of 3)", 1);
579 mTask = std::bind(&Normalizer::euler34, ph::_1, ph::_2);
582 this->cook(
"Normalizing level set using TVD_RK3 (step 2 of 3)", 2);
586 mTask = std::bind(&Normalizer::euler13, ph::_1, ph::_2);
589 this->cook(
"Normalizing level set using TVD_RK3 (step 3 of 3)", 2);
593 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
597 mTracker.mLeafs->removeAuxBuffers();
602template<
typename Gr
idT,
typename InterruptT>
607LevelSetTracker<GridT, InterruptT>::
608Normalizer<SpatialScheme, TemporalScheme, MaskT>::
609cook(
const char* msg,
int swapBuffer)
611 mTracker.startInterrupter( msg );
613 const int grainSize = mTracker.getGrainSize();
614 const LeafRange range = mTracker.leafs().leafRange(grainSize);
616 grainSize>0 ? tbb::parallel_for(range, *
this) : (*this)(range);
618 mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize==0);
620 mTracker.endInterrupter();
623template<
typename Gr
idT,
typename InterruptT>
627template <
int Nominator,
int Denominator>
629LevelSetTracker<GridT, InterruptT>::
630Normalizer<SpatialScheme, TemporalScheme, MaskT>::
631eval(StencilT& stencil,
const ValueType* phi, ValueType* result,
Index n)
const
633 using GradientT =
typename math::ISGradientNormSqrd<SpatialScheme>;
634 static const ValueType alpha = ValueType(Nominator)/ValueType(Denominator);
635 static const ValueType beta = ValueType(1) - alpha;
637 const ValueType normSqGradPhi = GradientT::result(stencil);
638 const ValueType phi0 = stencil.getValue();
641 v = phi0 - mDt * v * (
math::Sqrt(normSqGradPhi) * mInvDx - 1.0f);
642 result[n] = Nominator ? alpha * phi[n] + beta * v : v;
645template<
typename Gr
idT,
typename InterruptT>
649template <
int Nominator,
int Denominator>
651LevelSetTracker<GridT,InterruptT>::
652Normalizer<SpatialScheme, TemporalScheme, MaskT>::
653euler(
const LeafRange& range,
Index phiBuffer,
Index resultBuffer)
655 mTracker.checkInterrupter();
657 StencilT stencil(mTracker.grid());
659 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
660 const ValueType* phi = leafIter.buffer(phiBuffer).data();
661 ValueType* result = leafIter.buffer(resultBuffer).data();
662 if (mMask ==
nullptr) {
663 for (
auto iter = leafIter->cbeginValueOn(); iter; ++iter) {
664 stencil.moveTo(iter);
665 this->eval<Nominator, Denominator>(stencil, phi, result, iter.pos());
667 }
else if (
const MaskLeafT* mask = mMask->probeLeaf(leafIter->origin())) {
668 const ValueType* phi0 = leafIter->buffer().data();
669 for (MaskIterT iter = mask->cbeginValueOn(); iter; ++iter) {
670 const Index i = iter.pos();
671 stencil.moveTo(iter.getCoord(), phi0[i]);
672 this->eval<Nominator, Denominator>(stencil, phi, result, i);
684#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
686#ifdef OPENVDB_INSTANTIATE_LEVELSETTRACKER
Efficient multi-threaded replacement of the background values in tree.
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Implementation of morphological dilation and erosion.
Defined various multi-threaded utility functions for trees.
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition Types.h:683
Definition Exceptions.h:65
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition LeafManager.h:85
GridType
List of types that are currently supported by NanoVDB.
Definition NanoVDB.h:317
TemporalIntegrationScheme
Temporal integration schemes.
Definition FiniteDifference.h:233
@ TVD_RK1
Definition FiniteDifference.h:235
@ UNKNOWN_TIS
Definition FiniteDifference.h:234
@ TVD_RK2
Definition FiniteDifference.h:236
@ TVD_RK3
Definition FiniteDifference.h:237
float Sqrt(float x)
Return the square root of a floating-point value.
Definition Math.h:761
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition FiniteDifference.h:164
@ FIRST_BIAS
Definition FiniteDifference.h:166
@ THIRD_BIAS
Definition FiniteDifference.h:168
@ WENO5_BIAS
Definition FiniteDifference.h:169
@ UNKNOWN_BIAS
Definition FiniteDifference.h:165
@ SECOND_BIAS
Definition FiniteDifference.h:167
@ HJWENO5_BIAS
Definition FiniteDifference.h:170
float RoundDown(float x)
Return x rounded down to the nearest integer.
Definition Math.h:803
Type Pow2(Type x)
Return x2.
Definition Math.h:548
State
Definition IndexIterator.h:40
bool wasInterrupted(T *i, int percent=-1)
Definition NullInterrupter.h:49
Index32 Index
Definition Types.h:54
@ GRID_LEVEL_SET
Definition Types.h:455
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
Definition Operators.h:129
static T value()
Definition Math.h:148
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
#define OPENVDB_INSTANTIATE_CLASS
Definition version.h.in:153