39 #ifndef OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
40 #define OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
44 #include <openvdb/math/FiniteDifference.h>
65 template <
typename VelGr
idT,
typename Interpolator = BoxSampler>
72 DiscreteField(
const VelGridT &vel): mAccessor(vel.tree()), mTransform(&vel.transform()) {}
83 Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz), result);
90 return mAccessor.getValue(ijk);
94 const typename VelGridT::ConstAccessor mAccessor;
105 template <
typename ScalarT =
float>
121 inline VectorType operator() (
const Vec3d& xyz, ScalarType time)
const;
126 return (*
this)(ijk.asVec3d(), time);
170 template<
typename GridT,
171 typename FieldT = EnrightField<typename GridT::ValueType>,
186 mTracker(grid, interrupt), mField(field),
188 mTemporalScheme(math::
TVD_RK2) {}
229 size_t advect(ScalarType time0, ScalarType time1);
242 LevelSetAdvect(
const LevelSetAdvect& other);
244 LevelSetAdvect(LevelSetAdvect& other, tbb::split);
246 virtual ~LevelSetAdvect() {
if (mIsMaster) this->clearField();};
249 size_t advect(ScalarType time0, ScalarType time1);
251 void operator()(
const RangeType& r)
const
253 if (mTask) mTask(const_cast<LevelSetAdvect*>(
this), r);
257 void operator()(
const RangeType& r)
259 if (mTask) mTask(
this, r);
263 void join(
const LevelSetAdvect& other) { mMaxAbsV =
math::Max(mMaxAbsV, other.mMaxAbsV); }
265 typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
266 LevelSetAdvection& mParent;
268 const ScalarType mMinAbsV;
272 const bool mIsMaster;
274 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
276 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
278 typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
280 void sampleXformedField(
const RangeType& r, ScalarType time0, ScalarType time1);
281 void sampleAlignedField(
const RangeType& r, ScalarType time0, ScalarType time1);
283 void euler1(
const RangeType& r, ScalarType dt,
Index resultBuffer);
286 void euler2(
const RangeType& r, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer);
289 template<math::BiasedGradientScheme SpatialScheme>
290 size_t advect1(ScalarType time0, ScalarType time1);
294 size_t advect2(ScalarType time0, ScalarType time1);
299 size_t advect3(ScalarType time0, ScalarType time1);
308 void operator=(
const LevelSetAdvection& other) {}
312 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
316 switch (mSpatialScheme) {
318 return this->advect1<math::FIRST_BIAS >(time0, time1);
320 return this->advect1<math::SECOND_BIAS >(time0, time1);
322 return this->advect1<math::THIRD_BIAS >(time0, time1);
324 return this->advect1<math::WENO5_BIAS >(time0, time1);
326 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
333 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
334 template<math::BiasedGradientScheme SpatialScheme>
338 switch (mTemporalScheme) {
340 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
342 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
344 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
351 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
355 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
358 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
359 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
360 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
361 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
362 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
363 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
364 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
365 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
372 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
377 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
379 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
380 return tmp.advect(time0, time1);
385 template <
typename ScalarT>
386 inline math::Vec3<ScalarT>
389 static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
390 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
391 const ScalarT tr = cos(ScalarT(time) * phase);
392 const ScalarT a = sin(ScalarT(2.0)*Py);
393 const ScalarT b = -sin(ScalarT(2.0)*Px);
394 const ScalarT c = sin(ScalarT(2.0)*Pz);
396 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
405 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
415 mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()),
421 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
425 LevelSetAdvection<GridT, FieldT, InterruptT>::
426 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
427 LevelSetAdvect(
const LevelSetAdvect& other):
428 mParent(other.mParent),
430 mMinAbsV(other.mMinAbsV),
431 mMaxAbsV(other.mMaxAbsV),
438 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
442 LevelSetAdvection<GridT, FieldT, InterruptT>::
443 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
444 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
445 mParent(other.mParent),
447 mMinAbsV(other.mMinAbsV),
448 mMaxAbsV(other.mMaxAbsV),
455 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
461 advect(ScalarType time0, ScalarType time1)
465 const bool isForward = time0 < time1;
466 while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
468 mParent.mTracker.leafs().rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
470 const ScalarType dt = this->sampleField(time0, time1);
474 switch(TemporalScheme) {
478 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
480 this->cook(PARALLEL_FOR, 1);
485 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
487 this->cook(PARALLEL_FOR, 1);
491 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), 1, 1);
493 this->cook(PARALLEL_FOR, 1);
498 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
500 this->cook(PARALLEL_FOR, 1);
504 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), 1, 2);
506 this->cook(PARALLEL_FOR, 2);
510 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), 1, 2);
512 this->cook(PARALLEL_FOR, 2);
515 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
519 time0 += isForward ? dt : -dt;
521 mParent.mTracker.leafs().removeAuxBuffers();
524 mParent.mTracker.track();
529 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
532 inline typename GridT::ValueType
533 LevelSetAdvection<GridT, FieldT, InterruptT>::
534 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
535 sampleField(ScalarType time0, ScalarType time1)
538 const size_t leafCount = mParent.mTracker.leafs().leafCount();
539 if (leafCount==0)
return ScalarType(0.0);
540 mVec =
new VectorType*[leafCount];
541 if (mParent.mField.transform() == mParent.mTracker.grid().transform()) {
542 mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0, time1);
544 mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0, time1);
546 this->cook(PARALLEL_REDUCE);
548 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
549 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
550 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
551 const ScalarType dt =
math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
555 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
559 LevelSetAdvection<GridT, FieldT, InterruptT>::
560 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
561 sampleXformedField(
const RangeType& range, ScalarType time0, ScalarType time1)
563 const bool isForward = time0 < time1;
564 typedef typename LeafType::ValueOnCIter VoxelIterT;
565 const MapT& map = *mMap;
566 mParent.mTracker.checkInterrupter();
567 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
568 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
569 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
571 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
572 const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
574 vec[m] = isForward ? V : -V;
580 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
584 LevelSetAdvection<GridT, FieldT, InterruptT>::
585 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
586 sampleAlignedField(
const RangeType& range, ScalarType time0, ScalarType time1)
588 const bool isForward = time0 < time1;
589 typedef typename LeafType::ValueOnCIter VoxelIterT;
590 mParent.mTracker.checkInterrupter();
591 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
592 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
593 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
595 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
596 const VectorType V = mParent.mField(iter.getCoord(), time0);
598 vec[m] = isForward ? V : -V;
604 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
608 LevelSetAdvection<GridT, FieldT, InterruptT>::
609 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
612 if (mVec == NULL)
return;
613 for (
size_t n=0, e=mParent.mTracker.leafs().leafCount(); n<e; ++n)
delete [] mVec[n];
618 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
622 LevelSetAdvection<GridT, FieldT, InterruptT>::
623 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
624 cook(ThreadingMode mode,
size_t swapBuffer)
626 mParent.mTracker.startInterrupter(
"Advecting level set");
628 if (mParent.mTracker.getGrainSize()==0) {
629 (*this)(mParent.mTracker.leafs().getRange());
630 }
else if (mode == PARALLEL_FOR) {
631 tbb::parallel_for(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
632 }
else if (mode == PARALLEL_REDUCE) {
633 tbb::parallel_reduce(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
635 throw std::runtime_error(
"Undefined threading mode");
638 mParent.mTracker.leafs().swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
640 mParent.mTracker.endInterrupter();
645 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
649 LevelSetAdvection<GridT, FieldT, InterruptT>::
650 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
651 euler1(
const RangeType& range, ScalarType dt,
Index resultBuffer)
653 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
654 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
655 typedef typename LeafType::ValueOnCIter VoxelIterT;
656 mParent.mTracker.checkInterrupter();
657 const MapT& map = *mMap;
658 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
659 Stencil stencil(mParent.mTracker.grid());
660 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
661 BufferType& result = leafs.getBuffer(n, resultBuffer);
662 const VectorType* vec = mVec[n];
664 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
665 stencil.moveTo(iter);
667 result.setValue(iter.pos(), *iter - dt * V.dot(G));
674 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
678 LevelSetAdvection<GridT, FieldT, InterruptT>::
679 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
680 euler2(
const RangeType& range, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer)
682 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
683 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
684 typedef typename LeafType::ValueOnCIter VoxelIterT;
685 mParent.mTracker.checkInterrupter();
686 const MapT& map = *mMap;
687 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
688 const ScalarType beta = ScalarType(1.0) - alpha;
689 Stencil stencil(mParent.mTracker.grid());
690 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
691 const BufferType& phi = leafs.getBuffer(n, phiBuffer);
692 BufferType& result = leafs.getBuffer(n, resultBuffer);
693 const VectorType* vec = mVec[n];
695 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
696 stencil.moveTo(iter);
698 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
707 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
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
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:244
Definition: FiniteDifference.h:264
Type Pow2(Type x)
Return .
Definition: Math.h:458
Definition: Exceptions.h:88
Definition: FiniteDifference.h:196
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:505
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
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 math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk, const Vec3< typename Accessor::ValueType > &V)
Definition: Operators.h:810
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:351
Definition: FiniteDifference.h:198
Definition: FiniteDifference.h:195
Vec3< double > Vec3d
Definition: Vec3.h:625
Definition: FiniteDifference.h:263
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:274
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:194
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:566