37 #ifndef OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
38 #define OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
41 #include <openvdb/tools/GridTransformer.h>
42 #include <openvdb/math/FiniteDifference.h>
63 template <
typename VelGr
idT,
typename Interpolator = BoxSampler>
72 mAccessor(vel.getConstAccessor()),
73 mTransform(vel.transform()) {}
81 inline VectorType operator() (
const Vec3d& xyz, ScalarType)
const;
86 return mAccessor.getValue(ijk);
90 const AccessorType mAccessor;
97 template <
typename ScalarT =
float>
113 inline VectorType operator() (
const Vec3d& xyz, ScalarType time)
const;
118 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) {}
223 size_t advect(ScalarType time0, ScalarType time1);
235 LevelSetAdvect(
const LevelSetAdvect& other);
237 LevelSetAdvect(LevelSetAdvect& other, tbb::split);
239 virtual ~LevelSetAdvect() {
if (mIsMaster) this->clearField();};
242 size_t advect(ScalarType time0, ScalarType time1);
244 void operator()(
const RangeType& r)
const
246 if (mTask) mTask(const_cast<LevelSetAdvect*>(
this), r);
250 void operator()(
const RangeType& r)
252 if (mTask) mTask(
this, r);
256 void join(
const LevelSetAdvect& other) { mMaxAbsV =
math::Max(mMaxAbsV, other.mMaxAbsV); }
258 typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
259 LevelSetAdvection& mParent;
261 const ScalarType mMinAbsV;
265 const bool mIsMaster;
267 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
269 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
271 typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
273 void sampleXformedField(
const RangeType& r, ScalarType t);
274 void sampleAlignedField(
const RangeType& r, ScalarType t);
276 void euler1(
const RangeType& r, ScalarType dt,
Index resultBuffer);
279 void euler2(
const RangeType& r, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer);
282 template<math::BiasedGradientScheme SpatialScheme>
283 size_t advect1(ScalarType time0, ScalarType time1);
287 size_t advect2(ScalarType time0, ScalarType time1);
292 size_t advect3(ScalarType time0, ScalarType time1);
301 void operator=(
const LevelSetAdvection& other) {}
305 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
309 switch (mSpatialScheme) {
311 return this->advect1<math::FIRST_BIAS >(time0, time1);
313 return this->advect1<math::SECOND_BIAS >(time0, time1);
315 return this->advect1<math::THIRD_BIAS >(time0, time1);
317 return this->advect1<math::WENO5_BIAS >(time0, time1);
319 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
326 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
327 template<math::BiasedGradientScheme SpatialScheme>
331 switch (mTemporalScheme) {
333 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
335 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
337 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
344 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
348 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
351 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
352 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
353 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
354 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
355 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
356 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
357 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
358 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
365 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
370 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
372 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
373 return tmp.advect(time0, time1);
376 template <
typename VelGr
idT,
typename Interpolator>
377 inline typename VelGridT::ValueType
380 const VectorType coord = mTransform.worldToIndex(xyz);
382 Interpolator::template sample<AccessorType>(mAccessor, coord, result);
386 template <
typename ScalarT>
390 static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
391 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
392 const ScalarT tr = cos(ScalarT(time) * phase);
393 const ScalarT a = sin(ScalarT(2.0)*Py);
394 const ScalarT b = -sin(ScalarT(2.0)*Px);
395 const ScalarT c = sin(ScalarT(2.0)*Pz);
397 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
406 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
416 mMap(*(parent.mTracker.grid().transform().template constMap<MapT>())),
422 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
426 LevelSetAdvection<GridT, FieldT, InterruptT>::
427 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
428 LevelSetAdvect(
const LevelSetAdvect& other):
429 mParent(other.mParent),
431 mMinAbsV(other.mMinAbsV),
432 mMaxAbsV(other.mMaxAbsV),
439 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
443 LevelSetAdvection<GridT, FieldT, InterruptT>::
444 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
445 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
446 mParent(other.mParent),
448 mMinAbsV(other.mMinAbsV),
449 mMaxAbsV(other.mMaxAbsV),
456 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
462 advect(ScalarType time0, ScalarType time1)
465 while (time1 > time0 && mParent.mTracker.checkInterrupter()) {
467 mParent.mTracker.mLeafs->rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
469 const ScalarType dt = this->sampleField(time0, time1);
472 switch(TemporalScheme) {
476 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
478 this->cook(PARALLEL_FOR, 1);
483 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
485 this->cook(PARALLEL_FOR, 1);
489 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), 1, 1);
491 this->cook(PARALLEL_FOR, 1);
496 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
498 this->cook(PARALLEL_FOR, 1);
502 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), 1, 2);
504 this->cook(PARALLEL_FOR, 2);
508 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), 1, 2);
510 this->cook(PARALLEL_FOR, 2);
513 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
518 mParent.mTracker.mLeafs->removeAuxBuffers();
522 mParent.mTracker.track();
527 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
530 inline typename GridT::ValueType
531 LevelSetAdvection<GridT, FieldT, InterruptT>::
532 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
533 sampleField(ScalarType time0, ScalarType time1)
536 if (mParent.mTracker.mLeafs->leafCount()==0)
return ScalarType(0.0);
537 mVec =
new VectorType*[mParent.mTracker.mLeafs->leafCount()];
538 if (mParent.mField.transform() == mParent.mTracker.mGrid.transform()) {
539 mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0);
541 mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0);
543 this->cook(PARALLEL_REDUCE);
545 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
546 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) : ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
547 return math::Min(time1-time0, ScalarType(CFL*mParent.mTracker.voxelSize()/
math::Sqrt(mMaxAbsV)));
550 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
554 LevelSetAdvection<GridT, FieldT, InterruptT>::
555 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
556 sampleXformedField(
const RangeType& range, ScalarType time)
558 typedef typename LeafType::ValueOnCIter VoxelIterT;
559 const MapT& map =
mMap;
560 mParent.mTracker.checkInterrupter();
561 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
562 const LeafType& leaf = mParent.mTracker.mLeafs->leaf(n);
563 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
565 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
566 const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time);
574 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
578 LevelSetAdvection<GridT, FieldT, InterruptT>::
579 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
580 sampleAlignedField(
const RangeType& range, ScalarType time)
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.mLeafs->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(), time);
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.mLeafs->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.mLeafs->getRange());
623 }
else if (mode == PARALLEL_FOR) {
624 tbb::parallel_for(mParent.mTracker.mLeafs->getRange(mParent.mTracker.getGrainSize()), *
this);
625 }
else if (mode == PARALLEL_REDUCE) {
626 tbb::parallel_reduce(mParent.mTracker.mLeafs->getRange(mParent.mTracker.getGrainSize()), *
this);
628 throw std::runtime_error(
"Undefined threading mode");
631 mParent.mTracker.mLeafs->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 Stencil stencil(mParent.mTracker.mGrid);
652 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
653 BufferType& result = mParent.mTracker.mLeafs->getBuffer(n, resultBuffer);
654 const VectorType* vec = mVec[n];
656 for (VoxelIterT iter = mParent.mTracker.mLeafs->leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
657 stencil.moveTo(iter);
659 result.setValue(iter.pos(), *iter - dt * V.dot(G));
666 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
670 LevelSetAdvection<GridT, FieldT, InterruptT>::
671 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
672 euler2(
const RangeType& range, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer)
674 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
675 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
676 typedef typename LeafType::ValueOnCIter VoxelIterT;
677 mParent.mTracker.checkInterrupter();
678 const MapT& map =
mMap;
679 const ScalarType beta = ScalarType(1.0) - alpha;
680 Stencil stencil(mParent.mTracker.mGrid);
681 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
682 const BufferType& phi = mParent.mTracker.mLeafs->getBuffer(n, phiBuffer);
683 BufferType& result = mParent.mTracker.mLeafs->getBuffer(n, resultBuffer);
684 const VectorType* vec = mVec[n];
686 for (VoxelIterT iter = mParent.mTracker.mLeafs->leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
687 stencil.moveTo(iter);
689 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
698 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED