42 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 43 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 45 #include <tbb/parallel_for.h> 46 #include <boost/bind.hpp> 47 #include <boost/function.hpp> 48 #include <openvdb/Types.h> 49 #include <openvdb/math/Math.h> 50 #include <openvdb/util/NullInterrupter.h> 96 template<
typename VelocityGridT =
Vec3fGrid,
97 bool StaggeredVelocity =
false,
112 , mInterrupter(interrupter)
113 , mIntegrator( Scheme::
SEMI )
114 , mLimiter( Scheme::
CLAMP )
119 e.
add(velGrid.background().length());
120 mMaxVelocity = e.
max();
141 switch (mIntegrator) {
204 template<
typename VolumeGr
idT>
207 if (!inGrid.hasUniformVoxels()) {
210 const double d = mMaxVelocity*
math::Abs(dt)/inGrid.voxelSize()[0];
232 template<
typename VolumeGridT,
233 typename VolumeSamplerT>
234 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
236 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
237 const double dt = timeStep/mSubSteps;
238 const int n = this->getMaxDistance(inGrid, dt);
240 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
242 for (
int step = 1; step < mSubSteps; ++step) {
243 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
245 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
246 outGrid.swap( tmpGrid );
279 template<
typename VolumeGridT,
281 typename VolumeSamplerT>
282 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
const MaskGridT& mask,
double timeStep)
284 if (inGrid.transform() != mask.transform()) {
286 "resampling either of the two grids into the index space of the other.");
288 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
289 const double dt = timeStep/mSubSteps;
290 const int n = this->getMaxDistance(inGrid, dt);
292 outGrid->topologyIntersection( mask );
294 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
295 outGrid->topologyUnion( inGrid );
297 for (
int step = 1; step < mSubSteps; ++step) {
298 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
300 tmpGrid->topologyIntersection( mask );
302 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
303 tmpGrid->topologyUnion( inGrid );
304 outGrid.swap( tmpGrid );
314 void start(
const char* str)
const 316 if (mInterrupter) mInterrupter->start(str);
320 if (mInterrupter) mInterrupter->end();
322 bool interrupt()
const 325 tbb::task::self().cancel_group_execution();
331 template<
typename VolumeGr
idT,
typename VolumeSamplerT>
332 void cook(VolumeGridT& outGrid,
const VolumeGridT& inGrid,
double dt)
334 outGrid.tree().voxelizeActiveTiles();
335 switch (mIntegrator) {
337 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
338 adv.cook(outGrid, dt);
342 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
343 adv.cook(outGrid, dt);
347 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
348 adv.cook(outGrid, dt);
352 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
353 adv.cook(outGrid, dt);
357 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
358 adv.cook(outGrid, dt);
362 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
363 adv.cook(outGrid, dt);
373 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
376 const VelocityGridT& mVelGrid;
378 InterrupterType* mInterrupter;
386 template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
387 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
388 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
390 typedef typename VolumeGridT::TreeType TreeT;
391 typedef typename VolumeGridT::ConstAccessor AccT;
392 typedef typename TreeT::ValueType ValueT;
394 typedef typename LeafManagerT::LeafNodeType LeafNodeT;
395 typedef typename LeafManagerT::LeafRange LeafRangeT;
397 typedef typename VelocityIntegratorT::ElementType RealT;
398 typedef typename TreeT::LeafNodeType::ValueOnIter VoxelIterT;
404 , mVelocityInt(parent.mVelGrid)
408 inline void cook(
const LeafRangeT& range)
410 if (mParent->mGrainSize > 0) {
411 tbb::parallel_for(range, *
this);
416 void operator()(
const LeafRangeT& range)
const 419 mTask(const_cast<Advect*>(
this), range);
421 void cook(VolumeGridT& outGrid,
double time_step)
423 mParent->start(
"Advecting volume");
424 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
425 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
427 const RealT dt =
static_cast<RealT
>(-time_step);
429 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, *mInGrid);
431 mTask = boost::bind(&Advect::rk, _1, _2,-dt, 1, outGrid);
433 mTask = boost::bind(&Advect::mac, _1, _2);
436 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, *mInGrid);
438 mTask = boost::bind(&Advect::rk, _1, _2,-dt, 1, outGrid);
440 mTask = boost::bind(&Advect::bfecc, _1, _2);
442 mTask = boost::bind(&Advect::rk, _1, _2, dt, 1, outGrid);
444 manager.swapLeafBuffer(1);
446 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, *mInGrid);
450 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
452 mTask = boost::bind(&Advect::limiter, _1, _2, dt);
458 void mac(
const LeafRangeT& range)
const 460 if (mParent->interrupt())
return;
462 AccT acc = mInGrid->getAccessor();
463 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
464 ValueT* out0 = leafIter.buffer( 0 ).data();
465 const ValueT* out1 = leafIter.buffer( 1 ).data();
466 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
468 const ValueT* in0 = leaf->buffer().data();
469 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
470 const Index i = voxelIter.pos();
471 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
474 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
475 const Index i = voxelIter.pos();
476 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
482 void bfecc(
const LeafRangeT& range)
const 484 if (mParent->interrupt())
return;
486 AccT acc = mInGrid->getAccessor();
487 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
488 ValueT* out0 = leafIter.buffer( 0 ).data();
489 const ValueT* out1 = leafIter.buffer( 1 ).data();
490 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
492 const ValueT* in0 = leaf->buffer().data();
493 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
494 const Index i = voxelIter.pos();
495 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
498 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
499 const Index i = voxelIter.pos();
500 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
506 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT& grid)
const 508 if (mParent->interrupt())
return;
510 AccT acc = grid.getAccessor();
511 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
512 ValueT* phi = leafIter.buffer( n ).data();
513 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
514 ValueT& value = phi[voxelIter.pos()];
516 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
517 value = SamplerT::sample(acc, xform.
worldToIndex(wPos));
521 void limiter(
const LeafRangeT& range, RealT dt)
const 523 if (mParent->interrupt())
return;
524 const bool doLimiter = mParent->isLimiterOn();
526 ValueT data[2][2][2], vMin, vMax;
528 AccT acc = mInGrid->getAccessor();
529 const ValueT backg = mInGrid->background();
530 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
531 ValueT* phi = leafIter.buffer( 0 ).data();
532 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
533 ValueT& value = phi[voxelIter.pos()];
536 assert(OrderRK == 1);
538 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
540 Coord ijk = Coord::floor( iPos );
541 BoxSampler::getValues(data, acc, ijk);
545 }
else if (value < vMin || value > vMax ) {
546 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
547 value = BoxSampler::trilinearInterpolation( data, iPos );
553 leafIter->setValueOff( voxelIter.pos() );
559 typename boost::function<void (Advect*, const LeafRangeT&)> mTask;
560 const VolumeGridT* mInGrid;
561 const VelocityIntegratorT mVelocityInt;
569 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
Index32 Index
Definition: Types.h:58
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
This class computes the minimum and maximum values of a population of floating-point values...
Definition: Stats.h:55
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:293
Delta for small floating-point offsets.
Definition: Math.h:132
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:735
Definition: Exceptions.h:88
math::Vec3< Real > Vec3R
Definition: Types.h:76
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean, variance, etc.) of grid values.
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:114
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:370
Defined various multi-threaded utility functions for trees.
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Vec3SGrid Vec3fGrid
Definition: openvdb.h:77
Definition: Exceptions.h:86
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:39
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
Vec3< double > Vec3d
Definition: Vec3.h:643
void add(double val)
Add a single sample.
Definition: Stats.h:69
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:246
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
double max() const
Return the maximum value.
Definition: Stats.h:91