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;
101 template <
typename ScalarT =
float>
117 inline VectorType operator() (
const Vec3d& xyz, ScalarType time)
const;
122 return (*
this)(ijk.
asVec3d(), time);
166 template<
typename GridT,
167 typename FieldT = EnrightField<typename GridT::ValueType>,
182 mTracker(grid, interrupt), mField(field),
184 mTemporalScheme(math::
TVD_RK2) {}
225 size_t advect(ScalarType time0, ScalarType time1);
238 LevelSetAdvect(
const LevelSetAdvect& other);
240 LevelSetAdvect(LevelSetAdvect& other, tbb::split);
242 virtual ~LevelSetAdvect() {
if (mIsMaster) this->clearField();};
245 size_t advect(ScalarType time0, ScalarType time1);
247 void operator()(
const RangeType& r)
const
249 if (mTask) mTask(const_cast<LevelSetAdvect*>(
this), r);
253 void operator()(
const RangeType& r)
255 if (mTask) mTask(
this, r);
259 void join(
const LevelSetAdvect& other) { mMaxAbsV =
math::Max(mMaxAbsV, other.mMaxAbsV); }
261 typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
262 LevelSetAdvection& mParent;
264 const ScalarType mMinAbsV;
268 const bool mIsMaster;
270 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
272 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
274 typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
276 void sampleXformedField(
const RangeType& r, ScalarType time0, ScalarType time1);
277 void sampleAlignedField(
const RangeType& r, ScalarType time0, ScalarType time1);
279 void euler1(
const RangeType& r, ScalarType dt,
Index resultBuffer);
282 void euler2(
const RangeType& r, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer);
285 template<math::BiasedGradientScheme SpatialScheme>
286 size_t advect1(ScalarType time0, ScalarType time1);
290 size_t advect2(ScalarType time0, ScalarType time1);
295 size_t advect3(ScalarType time0, ScalarType time1);
304 void operator=(
const LevelSetAdvection& other) {}
308 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
312 switch (mSpatialScheme) {
314 return this->advect1<math::FIRST_BIAS >(time0, time1);
316 return this->advect1<math::SECOND_BIAS >(time0, time1);
318 return this->advect1<math::THIRD_BIAS >(time0, time1);
320 return this->advect1<math::WENO5_BIAS >(time0, time1);
322 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
329 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
330 template<math::BiasedGradientScheme SpatialScheme>
334 switch (mTemporalScheme) {
336 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
338 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
340 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
347 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
351 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
354 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
355 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
356 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
357 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
358 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
359 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
360 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
361 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
368 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
373 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
375 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
376 return tmp.advect(time0, time1);
381 template <
typename ScalarT>
382 inline math::Vec3<ScalarT>
385 static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
386 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
387 const ScalarT tr = cos(ScalarT(time) * phase);
388 const ScalarT a = sin(ScalarT(2.0)*Py);
389 const ScalarT b = -sin(ScalarT(2.0)*Px);
390 const ScalarT c = sin(ScalarT(2.0)*Pz);
392 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
401 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
411 mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()),
417 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
421 LevelSetAdvection<GridT, FieldT, InterruptT>::
422 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
423 LevelSetAdvect(
const LevelSetAdvect& other):
424 mParent(other.mParent),
426 mMinAbsV(other.mMinAbsV),
427 mMaxAbsV(other.mMaxAbsV),
434 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
438 LevelSetAdvection<GridT, FieldT, InterruptT>::
439 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
440 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
441 mParent(other.mParent),
443 mMinAbsV(other.mMinAbsV),
444 mMaxAbsV(other.mMaxAbsV),
451 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
457 advect(ScalarType time0, ScalarType time1)
461 const bool isForward = time0 < time1;
462 while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
464 mParent.mTracker.leafs().rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
466 const ScalarType dt = this->sampleField(time0, time1);
469 switch(TemporalScheme) {
473 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
475 this->cook(PARALLEL_FOR, 1);
480 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
482 this->cook(PARALLEL_FOR, 1);
486 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), 1, 1);
488 this->cook(PARALLEL_FOR, 1);
493 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
495 this->cook(PARALLEL_FOR, 1);
499 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), 1, 2);
501 this->cook(PARALLEL_FOR, 2);
505 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), 1, 2);
507 this->cook(PARALLEL_FOR, 2);
510 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
512 time0 += isForward ? dt : -dt;
514 mParent.mTracker.leafs().removeAuxBuffers();
517 mParent.mTracker.track();
522 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
525 inline typename GridT::ValueType
526 LevelSetAdvection<GridT, FieldT, InterruptT>::
527 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
528 sampleField(ScalarType time0, ScalarType time1)
531 const size_t leafCount = mParent.mTracker.leafs().leafCount();
532 if (leafCount==0)
return ScalarType(0.0);
533 mVec =
new VectorType*[leafCount];
534 if (mParent.mField.transform() == mParent.mTracker.grid().transform()) {
535 mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0, time1);
537 mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0, time1);
539 this->cook(PARALLEL_REDUCE);
541 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
542 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
543 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
544 const ScalarType dt =
math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
548 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
552 LevelSetAdvection<GridT, FieldT, InterruptT>::
553 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
554 sampleXformedField(
const RangeType& range, ScalarType time0, ScalarType time1)
556 const bool isForward = time0 < time1;
557 typedef typename LeafType::ValueOnCIter VoxelIterT;
558 const MapT& map = *
mMap;
559 mParent.mTracker.checkInterrupter();
560 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
561 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
562 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
564 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
565 const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
567 vec[m] = isForward ? V : -V;
573 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
577 LevelSetAdvection<GridT, FieldT, InterruptT>::
578 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
579 sampleAlignedField(
const RangeType& range, ScalarType time0, ScalarType time1)
581 const bool isForward = time0 < time1;
582 typedef typename LeafType::ValueOnCIter VoxelIterT;
583 mParent.mTracker.checkInterrupter();
584 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
585 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
586 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
588 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
589 const VectorType V = mParent.mField(iter.getCoord(), time0);
591 vec[m] = isForward ? V : -V;
597 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
601 LevelSetAdvection<GridT, FieldT, InterruptT>::
602 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
605 if (mVec == NULL)
return;
606 for (
size_t n=0, e=mParent.mTracker.leafs().leafCount(); n<e; ++n)
delete [] mVec[n];
611 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
615 LevelSetAdvection<GridT, FieldT, InterruptT>::
616 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
617 cook(ThreadingMode mode,
size_t swapBuffer)
619 mParent.mTracker.startInterrupter(
"Advecting level set");
621 if (mParent.mTracker.getGrainSize()==0) {
622 (*this)(mParent.mTracker.leafs().getRange());
623 }
else if (mode == PARALLEL_FOR) {
624 tbb::parallel_for(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
625 }
else if (mode == PARALLEL_REDUCE) {
626 tbb::parallel_reduce(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
628 throw std::runtime_error(
"Undefined threading mode");
631 mParent.mTracker.leafs().swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
633 mParent.mTracker.endInterrupter();
638 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
642 LevelSetAdvection<GridT, FieldT, InterruptT>::
643 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
644 euler1(
const RangeType& range, ScalarType dt,
Index resultBuffer)
646 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
647 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
648 typedef typename LeafType::ValueOnCIter VoxelIterT;
649 mParent.mTracker.checkInterrupter();
650 const MapT& map = *
mMap;
651 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
652 Stencil stencil(mParent.mTracker.grid());
653 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
654 BufferType& result = leafs.getBuffer(n, resultBuffer);
655 const VectorType* vec = mVec[n];
657 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
658 stencil.moveTo(iter);
660 result.setValue(iter.pos(), *iter - dt * V.dot(G));
667 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
671 LevelSetAdvection<GridT, FieldT, InterruptT>::
672 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
673 euler2(
const RangeType& range, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer)
675 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
676 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
677 typedef typename LeafType::ValueOnCIter VoxelIterT;
678 mParent.mTracker.checkInterrupter();
679 const MapT& map = *
mMap;
680 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
681 const ScalarType beta = ScalarType(1.0) - alpha;
682 Stencil stencil(mParent.mTracker.grid());
683 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
684 const BufferType& phi = leafs.getBuffer(n, phiBuffer);
685 BufferType& result = leafs.getBuffer(n, resultBuffer);
686 const VectorType* vec = mVec[n];
688 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
689 stencil.moveTo(iter);
691 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
700 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED