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>
74 mAccessor(vel.getConstAccessor()),
75 mTransform(vel.transform()) {}
83 inline VectorType operator() (
const Vec3d& xyz, ScalarType)
const;
88 return mAccessor.getValue(ijk);
92 const AccessorType mAccessor;
99 template <
typename ScalarT =
float>
115 inline VectorType operator() (
const Vec3d& xyz, ScalarType time)
const;
120 return (*
this)(ijk.
asVec3d(), time);
168 template<
typename GridT,
169 typename FieldT = EnrightField<typename GridT::ValueType>,
184 mTracker(grid, interrupt), mField(field),
186 mTemporalScheme(math::
TVD_RK2) {}
227 size_t advect(ScalarType time0, ScalarType time1);
240 LevelSetAdvect(
const LevelSetAdvect& other);
242 LevelSetAdvect(LevelSetAdvect& other, tbb::split);
244 virtual ~LevelSetAdvect() {
if (mIsMaster) this->clearField();};
247 size_t advect(ScalarType time0, ScalarType time1);
249 void operator()(
const RangeType& r)
const
251 if (mTask) mTask(const_cast<LevelSetAdvect*>(
this), r);
255 void operator()(
const RangeType& r)
257 if (mTask) mTask(
this, r);
261 void join(
const LevelSetAdvect& other) { mMaxAbsV =
math::Max(mMaxAbsV, other.mMaxAbsV); }
263 typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
264 LevelSetAdvection& mParent;
266 const ScalarType mMinAbsV;
270 const bool mIsMaster;
272 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
274 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
276 typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
278 void sampleXformedField(
const RangeType& r, ScalarType time0, ScalarType time1);
279 void sampleAlignedField(
const RangeType& r, ScalarType time0, ScalarType time1);
281 void euler1(
const RangeType& r, ScalarType dt,
Index resultBuffer);
284 void euler2(
const RangeType& r, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer);
287 template<math::BiasedGradientScheme SpatialScheme>
288 size_t advect1(ScalarType time0, ScalarType time1);
292 size_t advect2(ScalarType time0, ScalarType time1);
297 size_t advect3(ScalarType time0, ScalarType time1);
306 void operator=(
const LevelSetAdvection& other) {}
310 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
314 switch (mSpatialScheme) {
316 return this->advect1<math::FIRST_BIAS >(time0, time1);
318 return this->advect1<math::SECOND_BIAS >(time0, time1);
320 return this->advect1<math::THIRD_BIAS >(time0, time1);
322 return this->advect1<math::WENO5_BIAS >(time0, time1);
324 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
331 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
332 template<math::BiasedGradientScheme SpatialScheme>
336 switch (mTemporalScheme) {
338 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
340 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
342 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
349 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
353 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
356 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
357 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
358 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
359 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
360 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
361 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
362 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
363 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
370 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
375 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
377 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
378 return tmp.advect(time0, time1);
381 template <
typename VelGr
idT,
typename Interpolator>
382 inline typename VelGridT::ValueType
385 const VectorType coord = mTransform.worldToIndex(xyz);
387 Interpolator::template sample<AccessorType>(mAccessor, coord, result);
391 template <
typename ScalarT>
395 static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
396 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
397 const ScalarT tr = cos(ScalarT(time) * phase);
398 const ScalarT a = sin(ScalarT(2.0)*Py);
399 const ScalarT b = -sin(ScalarT(2.0)*Px);
400 const ScalarT c = sin(ScalarT(2.0)*Pz);
402 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
411 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
421 mMap(*(parent.mTracker.grid().transform().template constMap<MapT>())),
427 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
431 LevelSetAdvection<GridT, FieldT, InterruptT>::
432 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
433 LevelSetAdvect(
const LevelSetAdvect& other):
434 mParent(other.mParent),
436 mMinAbsV(other.mMinAbsV),
437 mMaxAbsV(other.mMaxAbsV),
444 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
448 LevelSetAdvection<GridT, FieldT, InterruptT>::
449 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
450 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
451 mParent(other.mParent),
453 mMinAbsV(other.mMinAbsV),
454 mMaxAbsV(other.mMaxAbsV),
461 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
467 advect(ScalarType time0, ScalarType time1)
471 const bool isForward = time0 < time1;
472 while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
474 mParent.mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
476 const ScalarType dt = this->sampleField(time0, time1);
479 switch(TemporalScheme) {
483 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
485 this->cook(PARALLEL_FOR, 1);
490 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
492 this->cook(PARALLEL_FOR, 1);
496 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), 1, 1);
498 this->cook(PARALLEL_FOR, 1);
503 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
505 this->cook(PARALLEL_FOR, 1);
509 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), 1, 2);
511 this->cook(PARALLEL_FOR, 2);
515 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), 1, 2);
517 this->cook(PARALLEL_FOR, 2);
520 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
522 time0 += isForward ? dt : -dt;
524 mParent.mTracker.mLeafs->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.mLeafs->leafCount();
542 if (leafCount==0)
return ScalarType(0.0);
543 mVec =
new VectorType*[leafCount];
544 if (mParent.mField.transform() == mParent.mTracker.mGrid.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 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
552 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
553 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
554 const ScalarType dt =
math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
558 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
562 LevelSetAdvection<GridT, FieldT, InterruptT>::
563 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
564 sampleXformedField(
const RangeType& range, ScalarType time0, ScalarType time1)
566 const bool isForward = time0 < time1;
567 typedef typename LeafType::ValueOnCIter VoxelIterT;
568 const MapT& map =
mMap;
569 mParent.mTracker.checkInterrupter();
570 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
571 const LeafType& leaf = mParent.mTracker.mLeafs->leaf(n);
572 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
574 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
575 const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
577 vec[m] = isForward ? V : -V;
583 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
587 LevelSetAdvection<GridT, FieldT, InterruptT>::
588 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
589 sampleAlignedField(
const RangeType& range, ScalarType time0, ScalarType time1)
591 const bool isForward = time0 < time1;
592 typedef typename LeafType::ValueOnCIter VoxelIterT;
593 mParent.mTracker.checkInterrupter();
594 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
595 const LeafType& leaf = mParent.mTracker.mLeafs->leaf(n);
596 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
598 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
599 const VectorType V = mParent.mField(iter.getCoord(), time0);
601 vec[m] = isForward ? V : -V;
607 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
611 LevelSetAdvection<GridT, FieldT, InterruptT>::
612 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
615 if (mVec == NULL)
return;
616 for (
size_t n=0, e=mParent.mTracker.mLeafs->leafCount(); n<e; ++n)
delete [] mVec[n];
621 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
625 LevelSetAdvection<GridT, FieldT, InterruptT>::
626 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
627 cook(ThreadingMode mode,
size_t swapBuffer)
629 mParent.mTracker.startInterrupter(
"Advecting level set");
631 if (mParent.mTracker.getGrainSize()==0) {
632 (*this)(mParent.mTracker.mLeafs->getRange());
633 }
else if (mode == PARALLEL_FOR) {
634 tbb::parallel_for(mParent.mTracker.mLeafs->getRange(mParent.mTracker.getGrainSize()), *
this);
635 }
else if (mode == PARALLEL_REDUCE) {
636 tbb::parallel_reduce(mParent.mTracker.mLeafs->getRange(mParent.mTracker.getGrainSize()), *
this);
638 throw std::runtime_error(
"Undefined threading mode");
641 mParent.mTracker.mLeafs->swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
643 mParent.mTracker.endInterrupter();
648 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
652 LevelSetAdvection<GridT, FieldT, InterruptT>::
653 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
654 euler1(
const RangeType& range, ScalarType dt,
Index resultBuffer)
656 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
657 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
658 typedef typename LeafType::ValueOnCIter VoxelIterT;
659 mParent.mTracker.checkInterrupter();
660 const MapT& map =
mMap;
661 Stencil stencil(mParent.mTracker.mGrid);
662 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
663 BufferType& result = mParent.mTracker.mLeafs->getBuffer(n, resultBuffer);
664 const VectorType* vec = mVec[n];
666 for (VoxelIterT iter = mParent.mTracker.mLeafs->leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
667 stencil.moveTo(iter);
669 result.setValue(iter.pos(), *iter - dt * V.dot(G));
676 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
680 LevelSetAdvection<GridT, FieldT, InterruptT>::
681 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
682 euler2(
const RangeType& range, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer)
684 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
685 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
686 typedef typename LeafType::ValueOnCIter VoxelIterT;
687 mParent.mTracker.checkInterrupter();
688 const MapT& map =
mMap;
689 const ScalarType beta = ScalarType(1.0) - alpha;
690 Stencil stencil(mParent.mTracker.mGrid);
691 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
692 const BufferType& phi = mParent.mTracker.mLeafs->getBuffer(n, phiBuffer);
693 BufferType& result = mParent.mTracker.mLeafs->getBuffer(n, resultBuffer);
694 const VectorType* vec = mVec[n];
696 for (VoxelIterT iter = mParent.mTracker.mLeafs->leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
697 stencil.moveTo(iter);
699 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
708 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED