13#ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
14#define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
19#include <openvdb/thread/Threading.h>
26#include <tbb/parallel_for.h>
72template<
typename VelocityGridT =
Vec3fGrid,
73 bool StaggeredVelocity =
false,
86 VolumeAdvection(
const VelocityGridT& velGrid, InterrupterType* interrupter =
nullptr)
88 , mInterrupter(interrupter)
89 , mIntegrator( Scheme::SEMI )
90 , mLimiter( Scheme::CLAMP )
95 e.
add(velGrid.background().length());
96 mMaxVelocity = e.
max();
109 mIntegrator == Scheme::BFECC) ? 2 : 1; }
117 switch (mIntegrator) {
118 case Scheme::SEMI:
return 1;
119 case Scheme::MID:
return 2;
120 case Scheme::RK3:
return 3;
121 case Scheme::RK4:
return 4;
122 case Scheme::BFECC:
return 2;
123 case Scheme::MAC:
return 2;
143 mLimiter != Scheme::NO_LIMITER; }
164 void setSubSteps(
int substeps) { mSubSteps = math::Max(1, substeps); }
180 template<
typename VolumeGr
idT>
183 if (!inGrid.hasUniformVoxels()) {
186 const double d = mMaxVelocity*math::Abs(dt)/inGrid.voxelSize()[0];
187 return static_cast<int>( math::RoundUp(d) );
208 template<
typename VolumeGridT,
209 typename VolumeSamplerT>
210 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
212 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
213 const double dt = timeStep/mSubSteps;
214 const int n = this->getMaxDistance(inGrid, dt);
216 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
217 for (
int step = 1; step < mSubSteps; ++step) {
218 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
220 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
221 outGrid.swap( tmpGrid );
254 template<
typename VolumeGridT,
256 typename VolumeSamplerT>
257 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
const MaskGridT& mask,
double timeStep)
259 if (inGrid.transform() != mask.transform()) {
261 "resampling either of the two grids into the index space of the other.");
263 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
264 const double dt = timeStep/mSubSteps;
265 const int n = this->getMaxDistance(inGrid, dt);
267 outGrid->topologyIntersection( mask );
269 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
270 outGrid->topologyUnion( inGrid );
272 for (
int step = 1; step < mSubSteps; ++step) {
273 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
275 tmpGrid->topologyIntersection( mask );
277 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
278 tmpGrid->topologyUnion( inGrid );
279 outGrid.swap( tmpGrid );
289 void start(
const char* str)
const
291 if (mInterrupter) mInterrupter->start(str);
295 if (mInterrupter) mInterrupter->end();
297 bool interrupt()
const
299 if (mInterrupter && util::wasInterrupted(mInterrupter)) {
300 thread::cancelGroupExecution();
306 template<
typename VolumeGr
idT,
typename VolumeSamplerT>
307 void cook(VolumeGridT& outGrid,
const VolumeGridT& inGrid,
double dt)
309 switch (mIntegrator) {
311 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
312 adv.cook(outGrid, dt);
316 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
317 adv.cook(outGrid, dt);
321 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
322 adv.cook(outGrid, dt);
326 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
327 adv.cook(outGrid, dt);
330 case Scheme::BFECC: {
331 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
332 adv.cook(outGrid, dt);
336 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
337 adv.cook(outGrid, dt);
341 OPENVDB_THROW(ValueError,
"Spatial difference scheme not supported!");
347 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
350 const VelocityGridT& mVelGrid;
352 InterrupterType* mInterrupter;
353 Scheme::SemiLagrangian mIntegrator;
354 Scheme::Limiter mLimiter;
360template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
361template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
362struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
364 using TreeT =
typename VolumeGridT::TreeType;
365 using AccT =
typename VolumeGridT::ConstAccessor;
366 using ValueT =
typename TreeT::ValueType;
367 using LeafManagerT =
typename tree::LeafManager<TreeT>;
368 using LeafNodeT =
typename LeafManagerT::LeafNodeType;
369 using LeafRangeT =
typename LeafManagerT::LeafRange;
370 using VelocityIntegratorT = VelocityIntegrator<VelocityGridT, StaggeredVelocity>;
371 using RealT =
typename VelocityIntegratorT::ElementType;
372 using VoxelIterT =
typename TreeT::LeafNodeType::ValueOnIter;
374 Advect(
const VolumeGridT& inGrid,
const VolumeAdvection& parent)
377 , mVelocityInt(parent.mVelGrid)
381 inline void cook(
const LeafRangeT& range)
383 if (mParent->mGrainSize > 0) {
384 tbb::parallel_for(range, *
this);
389 void operator()(
const LeafRangeT& range)
const
392 mTask(
const_cast<Advect*
>(
this), range);
394 void cook(VolumeGridT& outGrid,
double time_step)
396 namespace ph = std::placeholders;
398 mParent->start(
"Advecting volume");
399 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
400 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
401 const RealT dt =
static_cast<RealT
>(-time_step);
402 if (mParent->mIntegrator == Scheme::MAC) {
403 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
405 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
407 mTask = std::bind(&Advect::mac, ph::_1, ph::_2);
409 }
else if (mParent->mIntegrator == Scheme::BFECC) {
410 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
412 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
414 mTask = std::bind(&Advect::bfecc, ph::_1, ph::_2);
416 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 1, &outGrid);
418 manager.swapLeafBuffer(1);
420 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
424 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
426 mTask = std::bind(&Advect::limiter, ph::_1, ph::_2, dt);
432 void mac(
const LeafRangeT& range)
const
434 if (mParent->interrupt())
return;
435 assert( mParent->mIntegrator == Scheme::MAC );
436 AccT acc = mInGrid->getAccessor();
437 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
438 ValueT* out0 = leafIter.buffer( 0 ).data();
439 const ValueT* out1 = leafIter.buffer( 1 ).data();
440 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
441 if (leaf !=
nullptr) {
442 const ValueT* in0 = leaf->buffer().data();
443 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
444 const Index i = voxelIter.pos();
445 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
448 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
449 const Index i = voxelIter.pos();
450 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
456 void bfecc(
const LeafRangeT& range)
const
458 if (mParent->interrupt())
return;
459 assert( mParent->mIntegrator == Scheme::BFECC );
460 AccT acc = mInGrid->getAccessor();
461 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
462 ValueT* out0 = leafIter.buffer( 0 ).data();
463 const ValueT* out1 = leafIter.buffer( 1 ).data();
464 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
465 if (leaf !=
nullptr) {
466 const ValueT* in0 = leaf->buffer().data();
467 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
468 const Index i = voxelIter.pos();
469 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
472 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
473 const Index i = voxelIter.pos();
474 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
480 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT* grid)
const
482 if (mParent->interrupt())
return;
483 const math::Transform& xform = mInGrid->transform();
484 AccT acc = grid->getAccessor();
485 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
486 ValueT* phi = leafIter.buffer( n ).data();
487 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
488 ValueT& value = phi[voxelIter.pos()];
489 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
490 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
491 value = SamplerT::sample(acc, xform.worldToIndex(wPos));
495 void limiter(
const LeafRangeT& range, RealT dt)
const
497 if (mParent->interrupt())
return;
498 const bool doLimiter = mParent->isLimiterOn();
499 const bool doClamp = mParent->mLimiter == Scheme::CLAMP;
500 ValueT data[2][2][2], vMin, vMax;
501 const math::Transform& xform = mInGrid->transform();
502 AccT acc = mInGrid->getAccessor();
503 const ValueT backg = mInGrid->background();
504 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
505 ValueT* phi = leafIter.buffer( 0 ).data();
506 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
507 ValueT& value = phi[voxelIter.pos()];
510 assert(OrderRK == 1);
511 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
512 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
513 Vec3d iPos = xform.worldToIndex(wPos);
514 Coord ijk = Coord::floor( iPos );
515 BoxSampler::getValues(data, acc, ijk);
516 BoxSampler::extrema(data, vMin, vMax);
518 value = math::Clamp( value, vMin, vMax);
519 }
else if (value < vMin || value > vMax ) {
520 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
521 value = BoxSampler::trilinearInterpolation( data, iPos );
525 if (math::isApproxEqual(value, backg, math::Delta<ValueT>::value())) {
527 leafIter->setValueOff( voxelIter.pos() );
534 typename std::function<void (Advect*,
const LeafRangeT&)> mTask;
535 const VolumeGridT* mInGrid;
536 const VelocityIntegratorT mVelocityInt;
537 const VolumeAdvection* mParent;
546#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
548#ifdef OPENVDB_INSTANTIATE_VOLUMEADVECT
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Implementation of morphological dilation and erosion.
Defined various multi-threaded utility functions for trees.
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean,...
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
Container class that associates a tree with a transform and metadata.
Definition Grid.h:571
Definition Exceptions.h:63
This class computes the minimum and maximum values of a population of floating-point values.
Definition Stats.h:89
void add(double val)
Add a single sample.
Definition Stats.h:102
double max() const
Return the maximum value.
Definition Stats.h:124
Vec3< double > Vec3d
Definition Vec3.h:664
Index32 Index
Definition Types.h:54
Grid< FloatTree > FloatGrid
Definition openvdb.h:75
Grid< Vec3STree > Vec3SGrid
Definition openvdb.h:81
math::Vec3< Real > Vec3R
Definition Types.h:72
Grid< DoubleTree > DoubleGrid
Definition openvdb.h:74
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Base class for interrupters.
Definition NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
#define OPENVDB_INSTANTIATE_CLASS
Definition version.h.in:153
#define OPENVDB_INSTANTIATE
Definition version.h.in:152