38 #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED 39 #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED 43 #include <openvdb/math/FiniteDifference.h> 72 template<
typename GridT,
73 typename InterruptT = util::NullInterrupter>
86 LevelSetMorphing(GridT& sourceGrid,
const GridT& targetGrid, InterruptT* interrupt =
nullptr)
87 : mTracker(sourceGrid, interrupt)
88 , mTarget(&targetGrid)
91 , mTemporalScheme(math::
TVD_RK2)
101 void setTarget(
const GridT& targetGrid) { mTarget = &targetGrid; }
119 return mTracker.getSpatialScheme();
124 mTracker.setSpatialScheme(scheme);
129 return mTracker.getTemporalScheme();
134 mTracker.setTemporalScheme(scheme);
166 mDeltaMask = max-
min;
188 template<math::BiasedGradientScheme SpatialScheme>
201 const GridT *mTarget, *mMask;
215 Morph(
const Morph& other);
217 Morph(Morph& other, tbb::split);
224 void operator()(
const LeafRange& r)
const 226 if (mTask) mTask(const_cast<Morph*>(
this), r);
232 if (mTask) mTask(
this, r);
236 void join(
const Morph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
239 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
241 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
250 template <
int Nominator,
int Denominator>
253 inline void euler12(
const LeafRange& r,
ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
254 inline void euler34(
const LeafRange& r,
ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
255 inline void euler13(
const LeafRange& r,
ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
257 using FuncType =
typename std::function<void (Morph*, const LeafRange&)>;
266 template<
typename Gr
idT,
typename InterruptT>
270 switch (mSpatialScheme) {
272 return this->advect1<math::FIRST_BIAS >(time0, time1);
280 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
291 template<
typename Gr
idT,
typename InterruptT>
292 template<math::BiasedGradientScheme SpatialScheme>
296 switch (mTemporalScheme) {
298 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
300 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
302 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
310 template<
typename Gr
idT,
typename InterruptT>
317 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
318 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
319 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
320 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
322 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
323 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
324 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
325 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
332 template<
typename Gr
idT,
typename InterruptT>
339 Morph<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
340 return tmp.advect(time0, time1);
346 template<
typename Gr
idT,
typename InterruptT>
355 , mMap(parent.mTracker.grid().transform().template constMap<MapT>().
get())
360 template<
typename Gr
idT,
typename InterruptT>
366 Morph(
const Morph& other)
367 : mParent(other.mParent)
368 , mMinAbsS(other.mMinAbsS)
369 , mMaxAbsS(other.mMaxAbsS)
375 template<
typename Gr
idT,
typename InterruptT>
381 Morph(Morph& other, tbb::split)
382 : mParent(other.mParent)
383 , mMinAbsS(other.mMinAbsS)
384 , mMaxAbsS(other.mMaxAbsS)
390 template<
typename Gr
idT,
typename InterruptT>
398 namespace ph = std::placeholders;
404 while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
405 mParent->mTracker.
leafs().rebuildAuxBuffers(auxBuffers);
407 const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
411 switch(TemporalScheme) {
415 mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, 2);
418 this->cook(PARALLEL_FOR, 1);
423 mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, 2);
426 this->cook(PARALLEL_FOR, 1);
430 mTask = std::bind(&Morph::euler12, ph::_1, ph::_2, dt);
433 this->cook(PARALLEL_FOR, 1);
438 mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, 3);
441 this->cook(PARALLEL_FOR, 1);
445 mTask = std::bind(&Morph::euler34, ph::_1, ph::_2, dt);
448 this->cook(PARALLEL_FOR, 2);
452 mTask = std::bind(&Morph::euler13, ph::_1, ph::_2, dt);
455 this->cook(PARALLEL_FOR, 2);
465 mParent->mTracker.leafs().removeAuxBuffers();
468 mParent->mTracker.track();
474 template<
typename Gr
idT,
typename InterruptT>
477 inline typename GridT::ValueType
482 namespace ph = std::placeholders;
485 const size_t leafCount = mParent->mTracker.
leafs().leafCount();
486 if (leafCount==0 || time0 >= time1)
return ValueType(0);
489 if (mParent->mTarget->transform() == xform &&
490 (mParent->mMask ==
nullptr || mParent->mMask->transform() == xform)) {
491 mTask = std::bind(&Morph::sampleAlignedSpeed, ph::_1, ph::_2, speedBuffer);
493 mTask = std::bind(&Morph::sampleXformedSpeed, ph::_1, ph::_2, speedBuffer);
495 this->cook(PARALLEL_REDUCE);
500 const ValueType dt =
math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
501 return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
504 template<
typename Gr
idT,
typename InterruptT>
512 using VoxelIterT =
typename LeafType::ValueOnCIter;
515 const MapT& map = *mMap;
516 mParent->mTracker.checkInterrupter();
518 typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
519 SamplerT target(targetAcc, mParent->mTarget->transform());
520 if (mParent->mMask ==
nullptr) {
521 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
522 ValueType* speed = leafIter.buffer(speedBuffer).data();
524 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
526 s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
533 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
534 const bool invMask = mParent->isMaskInverted();
535 typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
536 SamplerT mask(maskAcc, mParent->mMask->transform());
537 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
538 ValueType* speed = leafIter.buffer(speedBuffer).data();
540 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
541 const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());
544 s -= target.wsSample(xyz);
545 s *= invMask ? 1 - a : a;
554 template<
typename Gr
idT,
typename InterruptT>
562 using VoxelIterT =
typename LeafType::ValueOnCIter;
566 typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
568 if (mParent->mMask ==
nullptr) {
569 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
570 ValueType* speed = leafIter.buffer(speedBuffer).data();
572 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
574 s -= target.getValue(voxelIter.getCoord());
581 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
582 const bool invMask = mParent->isMaskInverted();
583 typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
584 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
585 ValueType* speed = leafIter.buffer(speedBuffer).data();
587 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
588 const Coord ijk = voxelIter.getCoord();
591 s -= target.getValue(ijk);
592 s *= invMask ? 1 - a : a;
601 template<
typename Gr
idT,
typename InterruptT>
607 cook(ThreadingMode mode,
size_t swapBuffer)
611 const int grainSize = mParent->mTracker.getGrainSize();
612 const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
614 if (mParent->mTracker.getGrainSize()==0) {
616 }
else if (mode == PARALLEL_FOR) {
617 tbb::parallel_for(range, *
this);
618 }
else if (mode == PARALLEL_REDUCE) {
619 tbb::parallel_reduce(range, *
this);
622 <<
" or " <<
int(PARALLEL_REDUCE) <<
", got " <<
int(mode));
625 mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
627 mParent->mTracker.endInterrupter();
630 template<
typename Gr
idT,
typename InterruptT>
633 template <
int Nominator,
int Denominator>
641 using StencilT =
typename SchemeT::template ISStencil<GridType>::StencilType;
642 using VoxelIterT =
typename LeafType::ValueOnCIter;
648 mParent->mTracker.checkInterrupter();
649 const MapT& map = *mMap;
650 StencilT stencil(mParent->mTracker.grid());
652 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
653 const ValueType* speed = leafIter.buffer(speedBuffer).data();
655 const ValueType* phi = leafIter.buffer(phiBuffer).data();
656 ValueType* result = leafIter.buffer(resultBuffer).data();
657 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
658 const Index n = voxelIter.pos();
660 stencil.moveTo(voxelIter);
661 const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
662 result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
671 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
Definition: Operators.h:861
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:320
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:610
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:395
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:549
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:265
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:51
Definition: FiniteDifference.h:195
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: FiniteDifference.h:198
Definition: FiniteDifference.h:262
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x².
Definition: Math.h:256
Definition: FiniteDifference.h:193
Definition: Operators.h:153
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
Definition: Exceptions.h:92
Definition: Exceptions.h:40
Definition: FiniteDifference.h:194
Definition: FiniteDifference.h:263
Index32 Index
Definition: Types.h:61
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:715
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Definition: FiniteDifference.h:196
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:358
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:308
Definition: FiniteDifference.h:264