40 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 41 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 43 #include <tbb/parallel_for.h> 44 #include <boost/bind.hpp> 45 #include <boost/function.hpp> 46 #include <openvdb/Types.h> 47 #include <openvdb/math/Math.h> 48 #include <openvdb/util/NullInterrupter.h> 97 template<
typename VelocityGridT =
Vec3fGrid,
98 bool StaggeredVelocity =
false,
113 , mInterrupter(interrupter)
114 , mIntegrator( Scheme::
SEMI )
115 , mLimiter( Scheme::
CLAMP )
120 e.
add(velGrid.background().length());
121 mMaxVelocity = e.
max();
142 switch (mIntegrator) {
205 template<
typename VolumeGr
idT>
208 if (!inGrid.hasUniformVoxels()) {
211 const double d = mMaxVelocity*
math::Abs(dt)/inGrid.voxelSize()[0];
233 template<
typename VolumeGridT,
234 typename VolumeSamplerT>
235 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
237 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
238 const double dt = timeStep/mSubSteps;
239 const int n = this->getMaxDistance(inGrid, dt);
241 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 switch (mIntegrator) {
336 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
337 adv.cook(outGrid, dt);
341 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
342 adv.cook(outGrid, dt);
346 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
347 adv.cook(outGrid, dt);
351 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
352 adv.cook(outGrid, dt);
356 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
357 adv.cook(outGrid, dt);
361 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
362 adv.cook(outGrid, dt);
372 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
375 const VelocityGridT& mVelGrid;
377 InterrupterType* mInterrupter;
385 template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
386 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
387 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
389 typedef typename VolumeGridT::TreeType TreeT;
390 typedef typename VolumeGridT::ConstAccessor AccT;
391 typedef typename TreeT::ValueType ValueT;
393 typedef typename LeafManagerT::LeafNodeType LeafNodeT;
394 typedef typename LeafManagerT::LeafRange LeafRangeT;
396 typedef typename VelocityIntegratorT::ElementType RealT;
397 typedef typename TreeT::LeafNodeType::ValueOnIter VoxelIterT;
402 , mVelocityInt(parent.mVelGrid)
406 inline void cook(
const LeafRangeT& range)
408 if (mParent->mGrainSize > 0) {
409 tbb::parallel_for(range, *
this);
414 void operator()(
const LeafRangeT& range)
const 417 mTask(const_cast<Advect*>(
this), range);
419 void cook(VolumeGridT& outGrid,
double time_step)
421 mParent->start(
"Advecting volume");
422 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
423 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
424 const RealT dt =
static_cast<RealT
>(-time_step);
426 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, mInGrid);
428 mTask = boost::bind(&Advect::rk, _1, _2,-dt, 1, &outGrid);
430 mTask = boost::bind(&Advect::mac, _1, _2);
433 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, mInGrid);
435 mTask = boost::bind(&Advect::rk, _1, _2,-dt, 1, &outGrid);
437 mTask = boost::bind(&Advect::bfecc, _1, _2);
439 mTask = boost::bind(&Advect::rk, _1, _2, dt, 1, &outGrid);
441 manager.swapLeafBuffer(1);
443 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, mInGrid);
447 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
449 mTask = boost::bind(&Advect::limiter, _1, _2, dt);
455 void mac(
const LeafRangeT& range)
const 457 if (mParent->interrupt())
return;
459 AccT acc = mInGrid->getAccessor();
460 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
461 ValueT* out0 = leafIter.buffer( 0 ).data();
462 const ValueT* out1 = leafIter.buffer( 1 ).data();
463 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
465 const ValueT* in0 = leaf->buffer().data();
466 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
467 const Index i = voxelIter.pos();
468 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
471 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
472 const Index i = voxelIter.pos();
473 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
479 void bfecc(
const LeafRangeT& range)
const 481 if (mParent->interrupt())
return;
483 AccT acc = mInGrid->getAccessor();
484 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
485 ValueT* out0 = leafIter.buffer( 0 ).data();
486 const ValueT* out1 = leafIter.buffer( 1 ).data();
487 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
489 const ValueT* in0 = leaf->buffer().data();
490 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
491 const Index i = voxelIter.pos();
492 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
495 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
496 const Index i = voxelIter.pos();
497 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
503 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT* grid)
const 505 if (mParent->interrupt())
return;
507 AccT acc = grid->getAccessor();
508 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
509 ValueT* phi = leafIter.buffer( n ).data();
510 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
511 ValueT& value = phi[voxelIter.pos()];
513 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
514 value = SamplerT::sample(acc, xform.
worldToIndex(wPos));
518 void limiter(
const LeafRangeT& range, RealT dt)
const 520 if (mParent->interrupt())
return;
521 const bool doLimiter = mParent->isLimiterOn();
523 ValueT data[2][2][2], vMin, vMax;
525 AccT acc = mInGrid->getAccessor();
526 const ValueT backg = mInGrid->background();
527 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
528 ValueT* phi = leafIter.buffer( 0 ).data();
529 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
530 ValueT& value = phi[voxelIter.pos()];
533 assert(OrderRK == 1);
535 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
538 BoxSampler::getValues(data, acc, ijk);
542 }
else if (value < vMin || value > vMax ) {
543 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
544 value = BoxSampler::trilinearInterpolation( data, iPos );
550 leafIter->setValueOff( voxelIter.pos() );
557 typename boost::function<void (Advect*, const LeafRangeT&)> mTask;
558 const VolumeGridT* mInGrid;
559 const VelocityIntegratorT mVelocityInt;
567 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED double max() const
Return the maximum value.
Definition: Stats.h:141
Coord Abs(const Coord &xyz)
Definition: Coord.h:247
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:753
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
Delta for small floating-point offsets.
Definition: Math.h:132
#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:115
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:47
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:246
Defined various multi-threaded utility functions for trees.
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
Index32 Index
Definition: Types.h:58
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Vec3< double > Vec3d
Definition: Vec3.h:651
Definition: Exceptions.h:86
Implementation of morphological dilation and erosion.
static Coord floor(const Vec3< T > &xyz)
Return the largest integer coordinates that are not greater than xyz (node centered conversion)...
Definition: Coord.h:82
Definition: Exceptions.h:39
void add(double val)
Add a single sample.
Definition: Stats.h:119
This class computes the minimum and maximum values of a population of floating-point values...
Definition: Stats.h:105
Vec3SGrid Vec3fGrid
Definition: openvdb.h:80
Definition: Exceptions.h:88
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
math::Vec3< Real > Vec3R
Definition: Types.h:76
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76