40 #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED 41 #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED 45 #include <openvdb/math/FiniteDifference.h> 72 template<
typename GridT,
73 typename InterruptT = util::NullInterrupter>
87 const GridT& targetGrid,
88 InterruptT* interrupt = NULL)
89 : mTracker(sourceGrid, interrupt)
90 , mTarget(&targetGrid)
93 , mTemporalScheme(math::
TVD_RK2)
103 void setTarget(
const GridT& targetGrid) { mTarget = &targetGrid; }
121 return mTracker.getSpatialScheme();
126 mTracker.setSpatialScheme(scheme);
131 return mTracker.getTemporalScheme();
136 mTracker.setTemporalScheme(scheme);
151 ValueType
minMask()
const {
return mMinMask; }
155 ValueType
maxMask()
const {
return mDeltaMask + mMinMask; }
168 mDeltaMask = max-
min;
182 size_t advect(ValueType time0, ValueType time1);
190 template<math::BiasedGradientScheme SpatialScheme>
191 size_t advect1(ValueType time0, ValueType time1);
195 size_t advect2(ValueType time0, ValueType time1);
200 size_t advect3(ValueType time0, ValueType time1);
203 const GridT *mTarget, *mMask;
206 ValueType mMinMask, mDeltaMask;
217 Morph(
const Morph& other);
219 Morph(Morph& other, tbb::split);
224 size_t advect(ValueType time0, ValueType time1);
226 void operator()(
const LeafRange& r)
const 228 if (mTask) mTask(const_cast<Morph*>(
this), r);
232 void operator()(
const LeafRange& r)
234 if (mTask) mTask(
this, r);
238 void join(
const Morph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
241 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
243 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
246 typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1,
Index speedBuffer);
247 void sampleXformedSpeed(
const LeafRange& r,
Index speedBuffer);
248 void sampleAlignedSpeed(
const LeafRange& r,
Index speedBuffer);
252 template <
int Nominator,
int Denominator>
254 inline void euler01(
const LeafRange& r, ValueType t,
Index s) {this->euler<0,1>(r,t,0,1,s);}
255 inline void euler12(
const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
256 inline void euler34(
const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
257 inline void euler13(
const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
259 typedef typename boost::function<void (Morph*, const LeafRange&)> FuncType;
261 ValueType mMinAbsS, mMaxAbsS;
268 template<
typename Gr
idT,
typename InterruptT>
272 switch (mSpatialScheme) {
274 return this->advect1<math::FIRST_BIAS >(time0, time1);
282 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
289 template<
typename Gr
idT,
typename InterruptT>
290 template<math::BiasedGradientScheme SpatialScheme>
294 switch (mTemporalScheme) {
296 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
298 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
300 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
307 template<
typename Gr
idT,
typename InterruptT>
314 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
315 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
316 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
317 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
319 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
320 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
321 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
322 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
329 template<
typename Gr
idT,
typename InterruptT>
336 Morph<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
337 return tmp.advect(time0, time1);
343 template<
typename Gr
idT,
typename InterruptT>
352 , mMap(parent.mTracker.grid().transform().template constMap<MapT>().
get())
357 template<
typename Gr
idT,
typename InterruptT>
363 Morph(
const Morph& other)
364 : mParent(other.mParent)
365 , mMinAbsS(other.mMinAbsS)
366 , mMaxAbsS(other.mMaxAbsS)
372 template<
typename Gr
idT,
typename InterruptT>
378 Morph(Morph& other, tbb::split)
379 : mParent(other.mParent)
380 , mMinAbsS(other.mMinAbsS)
381 , mMaxAbsS(other.mMaxAbsS)
387 template<
typename Gr
idT,
typename InterruptT>
399 while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
402 const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
406 switch(TemporalScheme) {
410 mTask = boost::bind(&Morph::euler01, _1, _2, dt, 2);
413 this->cook(PARALLEL_FOR, 1);
418 mTask = boost::bind(&Morph::euler01, _1, _2, dt, 2);
421 this->cook(PARALLEL_FOR, 1);
425 mTask = boost::bind(&Morph::euler12, _1, _2, dt);
428 this->cook(PARALLEL_FOR, 1);
433 mTask = boost::bind(&Morph::euler01, _1, _2, dt, 3);
436 this->cook(PARALLEL_FOR, 1);
440 mTask = boost::bind(&Morph::euler34, _1, _2, dt);
443 this->cook(PARALLEL_FOR, 2);
447 mTask = boost::bind(&Morph::euler13, _1, _2, dt);
450 this->cook(PARALLEL_FOR, 2);
459 mParent->mTracker.leafs().removeAuxBuffers();
462 mParent->mTracker.track();
468 template<
typename Gr
idT,
typename InterruptT>
471 inline typename GridT::ValueType
478 if (leafCount==0 || time0 >= time1)
return ValueType(0);
481 if (mParent->mTarget->transform() == xform &&
482 (mParent->mMask == NULL || mParent->mMask->transform() == xform)) {
483 mTask = boost::bind(&Morph::sampleAlignedSpeed, _1, _2, speedBuffer);
485 mTask = boost::bind(&Morph::sampleXformedSpeed, _1, _2, speedBuffer);
487 this->cook(PARALLEL_REDUCE);
492 const ValueType dt =
math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
493 return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
496 template<
typename Gr
idT,
typename InterruptT>
504 typedef typename LeafType::ValueOnCIter VoxelIterT;
506 const MapT& map = *mMap;
507 mParent->mTracker.checkInterrupter();
509 typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
510 SamplerT target(targetAcc, mParent->mTarget->transform());
511 if (mParent->mMask == NULL) {
513 ValueType* speed = leafIter.buffer(speedBuffer).data();
515 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
517 s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
524 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
525 const bool invMask = mParent->isMaskInverted();
526 typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
527 SamplerT mask(maskAcc, mParent->mMask->transform());
529 ValueType* speed = leafIter.buffer(speedBuffer).data();
531 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
532 const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());
535 s -= target.wsSample(xyz);
536 s *= invMask ? 1 - a : a;
545 template<
typename Gr
idT,
typename InterruptT>
553 typedef typename LeafType::ValueOnCIter VoxelIterT;
556 typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
558 if (mParent->mMask == NULL) {
560 ValueType* speed = leafIter.buffer(speedBuffer).data();
562 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
564 s -= target.getValue(voxelIter.getCoord());
571 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
572 const bool invMask = mParent->isMaskInverted();
573 typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
575 ValueType* speed = leafIter.buffer(speedBuffer).data();
577 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
578 const Coord ijk = voxelIter.getCoord();
581 s -= target.getValue(ijk);
582 s *= invMask ? 1 - a : a;
591 template<
typename Gr
idT,
typename InterruptT>
597 cook(ThreadingMode mode,
size_t swapBuffer)
601 const int grainSize = mParent->mTracker.getGrainSize();
602 const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
604 if (mParent->mTracker.getGrainSize()==0) {
606 }
else if (mode == PARALLEL_FOR) {
607 tbb::parallel_for(range, *
this);
608 }
else if (mode == PARALLEL_REDUCE) {
609 tbb::parallel_reduce(range, *
this);
611 throw std::runtime_error(
"Undefined threading mode");
614 mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
616 mParent->mTracker.endInterrupter();
619 template<
typename Gr
idT,
typename InterruptT>
622 template <
int Nominator,
int Denominator>
630 typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
631 typedef typename LeafType::ValueOnCIter VoxelIterT;
637 mParent->mTracker.checkInterrupter();
638 const MapT& map = *mMap;
639 StencilT stencil(mParent->mTracker.grid());
642 const ValueType* speed = leafIter.buffer(speedBuffer).data();
644 const ValueType* phi = leafIter.buffer(phiBuffer).data();
645 ValueType* result = leafIter.buffer(resultBuffer).data();
646 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
647 const Index n = voxelIter.pos();
649 stencil.moveTo(voxelIter);
650 const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
651 result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
660 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
Coord Abs(const Coord &xyz)
Definition: Coord.h:247
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
Definition: LeafManager.h:132
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
void rebuildAuxBuffers(size_t auxBuffersPerLeaf, bool serial=false)
Change the number of auxiliary buffers.
Definition: LeafManager.h:314
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:727
Definition: Operators.h:860
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:622
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Definition: FiniteDifference.h:194
Definition: FiniteDifference.h:263
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:47
Definition: FiniteDifference.h:198
Definition: Operators.h:152
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:272
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
Definition: Exceptions.h:39
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Definition: FiniteDifference.h:265
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:335
Iterator begin() const
Definition: LeafManager.h:188
Definition: LeafManager.h:135
Definition: Exceptions.h:88
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:336
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Definition: FiniteDifference.h:264