39 #ifndef OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
40 #define OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
42 #include <openvdb/Platform.h>
45 #include <openvdb/math/FiniteDifference.h>
46 #include <boost/math/constants/constants.hpp>
67 template <
typename VelGr
idT,
typename Interpolator = BoxSampler>
74 DiscreteField(
const VelGridT &vel): mAccessor(vel.tree()), mTransform(&vel.transform()) {}
82 inline VectorType operator() (
const Vec3d& xyz, ScalarType)
const
84 VectorType result = zeroVal<VectorType>();
85 Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz), result);
90 inline VectorType operator() (
const Coord& ijk, ScalarType)
const
92 return mAccessor.getValue(ijk);
96 const typename VelGridT::ConstAccessor mAccessor;
107 template <
typename ScalarT =
float>
123 inline VectorType operator() (
const Vec3d& xyz, ScalarType time)
const;
126 inline VectorType operator() (
const Coord& ijk, ScalarType time)
const
128 return (*
this)(ijk.asVec3d(), time);
172 template<
typename GridT,
173 typename FieldT = EnrightField<typename GridT::ValueType>,
188 mTracker(grid, interrupt), mField(field),
190 mTemporalScheme(math::
TVD_RK2) {}
231 size_t advect(ScalarType time0, ScalarType time1);
244 LevelSetAdvect(
const LevelSetAdvect& other);
246 LevelSetAdvect(LevelSetAdvect& other, tbb::split);
248 virtual ~LevelSetAdvect() {
if (mIsMaster) this->clearField();};
251 size_t advect(ScalarType time0, ScalarType time1);
253 void operator()(
const RangeType& r)
const
255 if (mTask) mTask(const_cast<LevelSetAdvect*>(
this), r);
259 void operator()(
const RangeType& r)
261 if (mTask) mTask(
this, r);
265 void join(
const LevelSetAdvect& other) { mMaxAbsV =
math::Max(mMaxAbsV, other.mMaxAbsV); }
267 typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
268 LevelSetAdvection& mParent;
270 const ScalarType mMinAbsV;
274 const bool mIsMaster;
276 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
278 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
280 typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
282 void sampleXformedField(
const RangeType& r, ScalarType time0, ScalarType time1);
283 void sampleAlignedField(
const RangeType& r, ScalarType time0, ScalarType time1);
285 void euler1(
const RangeType& r, ScalarType dt,
Index resultBuffer);
288 void euler2(
const RangeType& r, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer);
291 template<math::BiasedGradientScheme SpatialScheme>
292 size_t advect1(ScalarType time0, ScalarType time1);
296 size_t advect2(ScalarType time0, ScalarType time1);
301 size_t advect3(ScalarType time0, ScalarType time1);
310 void operator=(
const LevelSetAdvection& other) {}
314 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
318 switch (mSpatialScheme) {
320 return this->advect1<math::FIRST_BIAS >(time0, time1);
322 return this->advect1<math::SECOND_BIAS >(time0, time1);
324 return this->advect1<math::THIRD_BIAS >(time0, time1);
326 return this->advect1<math::WENO5_BIAS >(time0, time1);
328 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
335 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
336 template<math::BiasedGradientScheme SpatialScheme>
340 switch (mTemporalScheme) {
342 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
344 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
346 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
353 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
357 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
360 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
361 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
362 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
363 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
364 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
365 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
366 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
367 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
374 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
379 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
381 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
382 return tmp.advect(time0, time1);
387 template <
typename ScalarT>
388 inline math::Vec3<ScalarT>
391 const ScalarT pi = boost::math::constants::pi<ScalarT>();
392 const ScalarT phase = pi / ScalarT(3.0);
393 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
394 const ScalarT tr = cos(ScalarT(time) * phase);
395 const ScalarT a = sin(ScalarT(2.0)*Py);
396 const ScalarT b = -sin(ScalarT(2.0)*Px);
397 const ScalarT c = sin(ScalarT(2.0)*Pz);
399 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
408 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
417 mMinAbsV(ScalarType(1e-6)),
418 mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()),
424 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
428 LevelSetAdvection<GridT, FieldT, InterruptT>::
429 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
430 LevelSetAdvect(
const LevelSetAdvect& other):
431 mParent(other.mParent),
433 mMinAbsV(other.mMinAbsV),
434 mMaxAbsV(other.mMaxAbsV),
441 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
445 LevelSetAdvection<GridT, FieldT, InterruptT>::
446 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
447 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
448 mParent(other.mParent),
450 mMinAbsV(other.mMinAbsV),
451 mMaxAbsV(other.mMaxAbsV),
458 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
464 advect(ScalarType time0, ScalarType time1)
468 const bool isForward = time0 < time1;
469 while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
471 mParent.mTracker.leafs().rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
473 const ScalarType dt = this->sampleField(time0, time1);
477 switch(TemporalScheme) {
481 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
483 this->cook(PARALLEL_FOR, 1);
488 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
490 this->cook(PARALLEL_FOR, 1);
494 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), 1, 1);
496 this->cook(PARALLEL_FOR, 1);
501 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
503 this->cook(PARALLEL_FOR, 1);
507 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), 1, 2);
509 this->cook(PARALLEL_FOR, 2);
513 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), 1, 2);
515 this->cook(PARALLEL_FOR, 2);
518 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
522 time0 += isForward ? dt : -dt;
524 mParent.mTracker.leafs().removeAuxBuffers();
527 mParent.mTracker.track();
532 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
535 inline typename GridT::ValueType
536 LevelSetAdvection<GridT, FieldT, InterruptT>::
537 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
538 sampleField(ScalarType time0, ScalarType time1)
541 const size_t leafCount = mParent.mTracker.leafs().leafCount();
542 if (leafCount==0)
return ScalarType(0.0);
543 mVec =
new VectorType*[leafCount];
544 if (mParent.mField.transform() == mParent.mTracker.grid().transform()) {
545 mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0, time1);
547 mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0, time1);
549 this->cook(PARALLEL_REDUCE);
551 #ifndef _MSC_VER // Visual C++ doesn't guarantee thread-safe initialization of local statics
554 const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
555 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
556 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
557 const ScalarType dt =
math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
561 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
565 LevelSetAdvection<GridT, FieldT, InterruptT>::
566 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
567 sampleXformedField(
const RangeType& range, ScalarType time0, ScalarType time1)
569 const bool isForward = time0 < time1;
570 typedef typename LeafType::ValueOnCIter VoxelIterT;
571 const MapT& map = *mMap;
572 mParent.mTracker.checkInterrupter();
573 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
574 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
575 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
577 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
578 const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
580 vec[m] = isForward ? V : -V;
586 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
590 LevelSetAdvection<GridT, FieldT, InterruptT>::
591 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
592 sampleAlignedField(
const RangeType& range, ScalarType time0, ScalarType time1)
594 const bool isForward = time0 < time1;
595 typedef typename LeafType::ValueOnCIter VoxelIterT;
596 mParent.mTracker.checkInterrupter();
597 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
598 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
599 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
601 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
602 const VectorType V = mParent.mField(iter.getCoord(), time0);
604 vec[m] = isForward ? V : -V;
610 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
614 LevelSetAdvection<GridT, FieldT, InterruptT>::
615 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
618 if (mVec == NULL)
return;
619 for (
size_t n=0, e=mParent.mTracker.leafs().leafCount(); n<e; ++n)
delete [] mVec[n];
624 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
628 LevelSetAdvection<GridT, FieldT, InterruptT>::
629 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
630 cook(ThreadingMode mode,
size_t swapBuffer)
632 mParent.mTracker.startInterrupter(
"Advecting level set");
634 if (mParent.mTracker.getGrainSize()==0) {
635 (*this)(mParent.mTracker.leafs().getRange());
636 }
else if (mode == PARALLEL_FOR) {
637 tbb::parallel_for(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
638 }
else if (mode == PARALLEL_REDUCE) {
639 tbb::parallel_reduce(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
641 throw std::runtime_error(
"Undefined threading mode");
644 mParent.mTracker.leafs().swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
646 mParent.mTracker.endInterrupter();
651 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
655 LevelSetAdvection<GridT, FieldT, InterruptT>::
656 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
657 euler1(
const RangeType& range, ScalarType dt,
Index resultBuffer)
659 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
660 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
661 typedef typename LeafType::ValueOnCIter VoxelIterT;
662 mParent.mTracker.checkInterrupter();
663 const MapT& map = *mMap;
664 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
665 Stencil stencil(mParent.mTracker.grid());
666 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
667 BufferType& result = leafs.getBuffer(n, resultBuffer);
668 const VectorType* vec = mVec[n];
670 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
671 stencil.moveTo(iter);
673 result.setValue(iter.pos(), *iter - dt * V.dot(G));
680 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
684 LevelSetAdvection<GridT, FieldT, InterruptT>::
685 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
686 euler2(
const RangeType& range, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer)
688 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
689 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
690 typedef typename LeafType::ValueOnCIter VoxelIterT;
691 mParent.mTracker.checkInterrupter();
692 const MapT& map = *mMap;
693 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
694 const ScalarType beta = ScalarType(1.0) - alpha;
695 Stencil stencil(mParent.mTracker.grid());
696 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
697 const BufferType& phi = leafs.getBuffer(n, phiBuffer);
698 BufferType& result = leafs.getBuffer(n, resultBuffer);
699 const VectorType* vec = mVec[n];
701 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
702 stencil.moveTo(iter);
704 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
713 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:265
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:284
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
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:831
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Type Pow2(Type x)
Return .
Definition: Math.h:498
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:606
Definition: FiniteDifference.h:196
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
Definition: Exceptions.h:39
Definition: FiniteDifference.h:195
Vec3< double > Vec3d
Definition: Vec3.h:629
Definition: FiniteDifference.h:264
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:545
Definition: Exceptions.h:88
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:391
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Definition: FiniteDifference.h:198
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Definition: FiniteDifference.h:194
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:314
Definition: FiniteDifference.h:263
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261