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;
260 LevelSetMorphing* mParent;
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>
311 LevelSetMorphing<GridT, InterruptT>::advect2(ValueType time0, ValueType time1)
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>
334 LevelSetMorphing<GridT, InterruptT>::advect3(ValueType time0, ValueType time1)
336 Morph<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
337 return tmp.advect(time0, time1);
343 template<
typename Gr
idT,
typename InterruptT>
347 LevelSetMorphing<GridT, InterruptT>::
348 Morph<MapT, SpatialScheme, TemporalScheme>::
349 Morph(LevelSetMorphing<GridT, InterruptT>& parent)
351 , mMinAbsS(ValueType(1e-6))
352 , mMap(parent.mTracker.grid().transform().template constMap<MapT>().get())
357 template<
typename Gr
idT,
typename InterruptT>
361 LevelSetMorphing<GridT, InterruptT>::
362 Morph<MapT, SpatialScheme, TemporalScheme>::
363 Morph(
const Morph& other)
364 : mParent(other.mParent)
365 , mMinAbsS(other.mMinAbsS)
366 , mMaxAbsS(other.mMaxAbsS)
372 template<
typename Gr
idT,
typename InterruptT>
376 LevelSetMorphing<GridT, InterruptT>::
377 Morph<MapT, SpatialScheme, TemporalScheme>::
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>
393 advect(ValueType time0, ValueType time1)
399 while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
400 mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
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);
453 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
459 mParent->mTracker.leafs().removeAuxBuffers();
462 mParent->mTracker.track();
468 template<
typename Gr
idT,
typename InterruptT>
471 inline typename GridT::ValueType
472 LevelSetMorphing<GridT, InterruptT>::
473 Morph<MapT, SpatialScheme, TemporalScheme>::
474 sampleSpeed(ValueType time0, ValueType time1,
Index speedBuffer)
477 const size_t leafCount = mParent->mTracker.leafs().leafCount();
478 if (leafCount==0 || time0 >= time1)
return ValueType(0);
480 const math::Transform& xform = mParent->mTracker.grid().transform();
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);
489 static const ValueType CFL = (TemporalScheme ==
math::TVD_RK1 ? ValueType(0.3) :
490 TemporalScheme == math::
TVD_RK2 ? ValueType(0.9) :
491 ValueType(1.0))/math::
Sqrt(ValueType(3.0));
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>
500 LevelSetMorphing<GridT, InterruptT>::
501 Morph<MapT, SpatialScheme, TemporalScheme>::
502 sampleXformedSpeed(
const LeafRange& range,
Index speedBuffer)
504 typedef typename LeafType::ValueOnCIter VoxelIterT;
505 typedef tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler> SamplerT;
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) {
512 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
513 ValueType* speed = leafIter.buffer(speedBuffer).data();
515 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
516 ValueType& s = speed[voxelIter.pos()];
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());
528 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
529 ValueType* speed = leafIter.buffer(speedBuffer).data();
531 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
532 const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());
534 ValueType& s = speed[voxelIter.pos()];
535 s -= target.wsSample(xyz);
536 s *= invMask ? 1 - a : a;
545 template<
typename Gr
idT,
typename InterruptT>
549 LevelSetMorphing<GridT, InterruptT>::
550 Morph<MapT, SpatialScheme, TemporalScheme>::
551 sampleAlignedSpeed(
const LeafRange& range,
Index speedBuffer)
553 typedef typename LeafType::ValueOnCIter VoxelIterT;
554 mParent->mTracker.checkInterrupter();
556 typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
558 if (mParent->mMask == NULL) {
559 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
560 ValueType* speed = leafIter.buffer(speedBuffer).data();
562 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
563 ValueType& s = speed[voxelIter.pos()];
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();
574 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
575 ValueType* speed = leafIter.buffer(speedBuffer).data();
577 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
578 const Coord ijk = voxelIter.getCoord();
580 ValueType& s = speed[voxelIter.pos()];
581 s -= target.getValue(ijk);
582 s *= invMask ? 1 - a : a;
591 template<
typename Gr
idT,
typename InterruptT>
595 LevelSetMorphing<GridT, InterruptT>::
596 Morph<MapT, SpatialScheme, TemporalScheme>::
597 cook(ThreadingMode mode,
size_t swapBuffer)
599 mParent->mTracker.startInterrupter(
"Morphing level set");
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>
624 LevelSetMorphing<GridT,InterruptT>::
625 Morph<MapT, SpatialScheme, TemporalScheme>::
626 euler(
const LeafRange& range, ValueType dt,
629 typedef math::BIAS_SCHEME<SpatialScheme> SchemeT;
630 typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
631 typedef typename LeafType::ValueOnCIter VoxelIterT;
632 typedef math::GradientNormSqrd<MapT, SpatialScheme> NumGrad;
634 static const ValueType Alpha = ValueType(Nominator)/ValueType(Denominator);
635 static const ValueType Beta = ValueType(1) - Alpha;
637 mParent->mTracker.checkInterrupter();
638 const MapT& map = *mMap;
639 StencilT stencil(mParent->mTracker.grid());
641 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
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
Index32 Index
Definition: Types.h:58
Definition: FiniteDifference.h:198
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:293
Definition: FiniteDifference.h:263
Definition: Exceptions.h:88
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
math::Vec3< Real > Vec3R
Definition: Types.h:76
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
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
#define OPENVDB_VERSION_NAME
Definition: version.h:43
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
Definition: Exceptions.h:39
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
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
Definition: FiniteDifference.h:194
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:272
Definition: FiniteDifference.h:265
Definition: FiniteDifference.h:264
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:709
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:622
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: LeafManager.h:131