38 #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED 39 #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED 70 template<
typename GridT,
71 typename InterruptT = util::NullInterrupter>
84 LevelSetMorphing(GridT& sourceGrid,
const GridT& targetGrid, InterruptT* interrupt =
nullptr)
85 : mTracker(sourceGrid, interrupt)
86 , mTarget(&targetGrid)
89 , mTemporalScheme(math::
TVD_RK2)
99 void setTarget(
const GridT& targetGrid) { mTarget = &targetGrid; }
117 return mTracker.getSpatialScheme();
122 mTracker.setSpatialScheme(scheme);
127 return mTracker.getTemporalScheme();
132 mTracker.setTemporalScheme(scheme);
147 ValueType
minMask()
const {
return mMinMask; }
151 ValueType
maxMask()
const {
return mDeltaMask + mMinMask; }
164 mDeltaMask = max-
min;
178 size_t advect(ValueType time0, ValueType time1);
186 template<math::BiasedGradientScheme SpatialScheme>
187 size_t advect1(ValueType time0, ValueType time1);
191 size_t advect2(ValueType time0, ValueType time1);
196 size_t advect3(ValueType time0, ValueType time1);
199 const GridT *mTarget, *mMask;
202 ValueType mMinMask, mDeltaMask;
213 Morph(
const Morph& other);
215 Morph(Morph& other, tbb::split);
220 size_t advect(ValueType time0, ValueType time1);
222 void operator()(
const LeafRange& r)
const 224 if (mTask) mTask(const_cast<Morph*>(
this), r);
228 void operator()(
const LeafRange& r)
230 if (mTask) mTask(
this, r);
234 void join(
const Morph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
237 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
239 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
242 typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1,
Index speedBuffer);
243 void sampleXformedSpeed(
const LeafRange& r,
Index speedBuffer);
244 void sampleAlignedSpeed(
const LeafRange& r,
Index speedBuffer);
248 template <
int Nominator,
int Denominator>
250 inline void euler01(
const LeafRange& r, ValueType t,
Index s) {this->euler<0,1>(r,t,0,1,s);}
251 inline void euler12(
const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
252 inline void euler34(
const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
253 inline void euler13(
const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
255 typedef typename boost::function<void (Morph*, const LeafRange&)> FuncType;
257 ValueType mMinAbsS, mMaxAbsS;
264 template<
typename Gr
idT,
typename InterruptT>
268 switch (mSpatialScheme) {
270 return this->advect1<math::FIRST_BIAS >(time0, time1);
278 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);
308 template<
typename Gr
idT,
typename InterruptT>
315 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
316 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
317 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
318 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
320 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
321 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
322 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
323 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
330 template<
typename Gr
idT,
typename InterruptT>
337 Morph<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
338 return tmp.advect(time0, time1);
344 template<
typename Gr
idT,
typename InterruptT>
353 , mMap(parent.mTracker.grid().transform().template constMap<MapT>().
get())
358 template<
typename Gr
idT,
typename InterruptT>
364 Morph(
const Morph& other)
365 : mParent(other.mParent)
366 , mMinAbsS(other.mMinAbsS)
367 , mMaxAbsS(other.mMaxAbsS)
373 template<
typename Gr
idT,
typename InterruptT>
379 Morph(Morph& other, tbb::split)
380 : mParent(other.mParent)
381 , mMinAbsS(other.mMinAbsS)
382 , mMaxAbsS(other.mMaxAbsS)
388 template<
typename Gr
idT,
typename InterruptT>
400 while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
403 const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
407 switch(TemporalScheme) {
411 mTask = boost::bind(&Morph::euler01, _1, _2, dt, 2);
414 this->cook(PARALLEL_FOR, 1);
419 mTask = boost::bind(&Morph::euler01, _1, _2, dt, 2);
422 this->cook(PARALLEL_FOR, 1);
426 mTask = boost::bind(&Morph::euler12, _1, _2, dt);
429 this->cook(PARALLEL_FOR, 1);
434 mTask = boost::bind(&Morph::euler01, _1, _2, dt, 3);
437 this->cook(PARALLEL_FOR, 1);
441 mTask = boost::bind(&Morph::euler34, _1, _2, dt);
444 this->cook(PARALLEL_FOR, 2);
448 mTask = boost::bind(&Morph::euler13, _1, _2, dt);
451 this->cook(PARALLEL_FOR, 2);
461 mParent->mTracker.leafs().removeAuxBuffers();
464 mParent->mTracker.track();
470 template<
typename Gr
idT,
typename InterruptT>
473 inline typename GridT::ValueType
480 if (leafCount==0 || time0 >= time1)
return ValueType(0);
483 if (mParent->mTarget->transform() == xform &&
484 (mParent->mMask ==
nullptr || mParent->mMask->transform() == xform)) {
485 mTask = boost::bind(&Morph::sampleAlignedSpeed, _1, _2, speedBuffer);
487 mTask = boost::bind(&Morph::sampleXformedSpeed, _1, _2, speedBuffer);
489 this->cook(PARALLEL_REDUCE);
494 const ValueType dt =
math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
495 return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
498 template<
typename Gr
idT,
typename InterruptT>
506 typedef typename LeafType::ValueOnCIter VoxelIterT;
508 const MapT& map = *mMap;
509 mParent->mTracker.checkInterrupter();
511 typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
512 SamplerT target(targetAcc, mParent->mTarget->transform());
513 if (mParent->mMask ==
nullptr) {
515 ValueType* speed = leafIter.buffer(speedBuffer).data();
517 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
519 s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
526 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
527 const bool invMask = mParent->isMaskInverted();
528 typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
529 SamplerT mask(maskAcc, mParent->mMask->transform());
531 ValueType* speed = leafIter.buffer(speedBuffer).data();
533 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
534 const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());
537 s -= target.wsSample(xyz);
538 s *= invMask ? 1 - a : a;
547 template<
typename Gr
idT,
typename InterruptT>
555 typedef typename LeafType::ValueOnCIter VoxelIterT;
558 typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
560 if (mParent->mMask ==
nullptr) {
562 ValueType* speed = leafIter.buffer(speedBuffer).data();
564 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
566 s -= target.getValue(voxelIter.getCoord());
573 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
574 const bool invMask = mParent->isMaskInverted();
575 typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
577 ValueType* speed = leafIter.buffer(speedBuffer).data();
579 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
580 const Coord ijk = voxelIter.getCoord();
583 s -= target.getValue(ijk);
584 s *= invMask ? 1 - a : a;
593 template<
typename Gr
idT,
typename InterruptT>
599 cook(ThreadingMode mode,
size_t swapBuffer)
603 const int grainSize = mParent->mTracker.getGrainSize();
604 const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
606 if (mParent->mTracker.getGrainSize()==0) {
608 }
else if (mode == PARALLEL_FOR) {
609 tbb::parallel_for(range, *
this);
610 }
else if (mode == PARALLEL_REDUCE) {
611 tbb::parallel_reduce(range, *
this);
613 throw std::runtime_error(
"Undefined threading mode");
616 mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
618 mParent->mTracker.endInterrupter();
621 template<
typename Gr
idT,
typename InterruptT>
624 template <
int Nominator,
int Denominator>
632 typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
633 typedef typename LeafType::ValueOnCIter VoxelIterT;
639 mParent->mTracker.checkInterrupter();
640 const MapT& map = *mMap;
641 StencilT stencil(mParent->mTracker.grid());
644 const ValueType* speed = leafIter.buffer(speedBuffer).data();
646 const ValueType* phi = leafIter.buffer(phiBuffer).data();
647 ValueType* result = leafIter.buffer(resultBuffer).data();
648 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
649 const Index n = voxelIter.pos();
651 stencil.moveTo(voxelIter);
652 const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
653 result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
662 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
Definition: LeafManager.h:132
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
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:272
Coord Abs(const Coord &xyz)
Definition: Coord.h:254
Definition: FiniteDifference.h:262
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:101
Definition: Operators.h:152
Definition: FiniteDifference.h:195
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:333
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:727
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Definition: FiniteDifference.h:264
Definition: Exceptions.h:39
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: FiniteDifference.h:194
Definition: FiniteDifference.h:263
Definition: FiniteDifference.h:196
void rebuildAuxBuffers(size_t auxBuffersPerLeaf, bool serial=false)
Change the number of auxiliary buffers.
Definition: LeafManager.h:312
Definition: Exceptions.h:92
Definition: FiniteDifference.h:265
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
Definition: FiniteDifference.h:198
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
Definition: FiniteDifference.h:193
Definition: LeafManager.h:135
Definition: Operators.h:860
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Index32 Index
Definition: Types.h:57
Iterator begin() const
Definition: LeafManager.h:186
Definition: FiniteDifference.h:197
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324