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 ScalarType
minMask()
const {
return mMinMask; }
155 ScalarType
maxMask()
const {
return mDeltaMask + mMinMask; }
168 mDeltaMask = max-
min;
182 size_t advect(ScalarType time0, ScalarType time1);
186 template<math::BiasedGradientScheme SpatialScheme>
187 size_t advect1(ScalarType time0, ScalarType time1);
191 size_t advect2(ScalarType time0, ScalarType time1);
196 size_t advect3(ScalarType time0, ScalarType time1);
199 const GridT *mTarget, *mMask;
202 ScalarType mMinMask, mDeltaMask;
215 LevelSetMorph(LevelSetMorphing<GridT, InterruptT>& parent);
217 LevelSetMorph(
const LevelSetMorph& other);
219 LevelSetMorph(LevelSetMorph& other, tbb::split);
221 virtual ~LevelSetMorph() {}
224 size_t advect(ScalarType time0, ScalarType time1);
226 void operator()(
const LeafRange& r)
const
228 if (mTask) mTask(const_cast<LevelSetMorph*>(
this), r);
232 void operator()(
const LeafRange& r)
234 if (mTask) mTask(
this, r);
235 else OPENVDB_THROW(ValueError,
"task is undefined - don\'t call this method directly");
238 void join(
const LevelSetMorph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
240 typedef typename boost::function<void (LevelSetMorph*, const LeafRange&)> FuncType;
241 LevelSetMorphing* mParent;
242 ScalarType mMinAbsS, mMaxAbsS;
247 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
249 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
252 typename GridT::ValueType sampleSpeed(ScalarType time0, ScalarType time1,
Index speedBuffer);
253 void sampleXformedSpeed(
const LeafRange& r,
Index speedBuffer);
254 void sampleAlignedSpeed(
const LeafRange& r,
Index speedBuffer);
257 void euler1(
const LeafRange& r, ScalarType dt,
Index resultBuffer,
Index speedBuffer);
261 void euler2(
const LeafRange& r, ScalarType dt, ScalarType alpha,
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(ScalarType time0, ScalarType 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(ScalarType time0, ScalarType time1)
336 LevelSetMorph<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
337 return tmp.advect(time0, time1);
343 template<
typename Gr
idT,
typename InterruptT>
347 LevelSetMorphing<GridT, InterruptT>::
348 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
349 LevelSetMorph(LevelSetMorphing<GridT, InterruptT>& parent)
351 , mMinAbsS(ScalarType(1e-6))
352 , mMap(parent.mTracker.grid().transform().template constMap<MapT>().get())
357 template<
typename Gr
idT,
typename InterruptT>
361 LevelSetMorphing<GridT, InterruptT>::
362 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
363 LevelSetMorph(
const LevelSetMorph& 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 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
378 LevelSetMorph(LevelSetMorph& other, tbb::split)
379 : mParent(other.mParent)
380 , mMinAbsS(other.mMinAbsS)
381 , mMaxAbsS(other.mMaxAbsS)
387 template<
typename Gr
idT,
typename InterruptT>
393 advect(ScalarType time0, ScalarType time1)
399 while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
400 mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
402 const ScalarType dt = this->sampleSpeed(time0, time1, auxBuffers);
406 switch(TemporalScheme) {
410 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 2);
412 this->cook(PARALLEL_FOR, 1);
417 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 2);
419 this->cook(PARALLEL_FOR, 1);
423 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.5),
426 this->cook(PARALLEL_FOR, 1);
431 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 3);
433 this->cook(PARALLEL_FOR, 1);
437 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.75),
440 this->cook(PARALLEL_FOR, 2);
444 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(1.0/3.0),
447 this->cook(PARALLEL_FOR, 2);
450 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
456 mParent->mTracker.leafs().removeAuxBuffers();
459 mParent->mTracker.track();
465 template<
typename Gr
idT,
typename InterruptT>
468 inline typename GridT::ValueType
469 LevelSetMorphing<GridT, InterruptT>::
470 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
471 sampleSpeed(ScalarType time0, ScalarType time1,
Index speedBuffer)
474 const size_t leafCount = mParent->mTracker.leafs().leafCount();
475 if (leafCount==0 || time0 >= time1)
return ScalarType(0);
477 const math::Transform& xform = mParent->mTracker.grid().transform();
478 if (mParent->mTarget->transform() == xform &&
479 (mParent->mMask == NULL || mParent->mMask->transform() == xform)) {
480 mTask = boost::bind(&LevelSetMorph::sampleAlignedSpeed, _1, _2, speedBuffer);
482 mTask = boost::bind(&LevelSetMorph::sampleXformedSpeed, _1, _2, speedBuffer);
484 this->cook(PARALLEL_REDUCE);
486 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
487 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
488 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
489 const ScalarType dt =
math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
490 return math::Min(dt, ScalarType(CFL*dx/mMaxAbsS));
493 template<
typename Gr
idT,
typename InterruptT>
497 LevelSetMorphing<GridT, InterruptT>::
498 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
499 sampleXformedSpeed(
const LeafRange& range,
Index speedBuffer)
501 typedef typename LeafType::ValueOnCIter VoxelIterT;
502 typedef tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler> SamplerT;
503 const MapT& map = *mMap;
504 mParent->mTracker.checkInterrupter();
506 typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
507 SamplerT target(targetAcc, mParent->mTarget->transform());
508 if (mParent->mMask == NULL) {
509 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
510 BufferType& speed = leafIter.buffer(speedBuffer);
511 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
512 ScalarType& s =
const_cast<ScalarType&
>(speed.getValue(voxelIter.pos()));
513 s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
518 const ScalarType
min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
519 const bool invMask = mParent->isMaskInverted();
520 typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
521 SamplerT mask(maskAcc, mParent->mMask->transform());
522 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
523 BufferType& source = leafIter.buffer(speedBuffer);
524 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
525 const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());
527 ScalarType& s =
const_cast<ScalarType&
>(source.getValue(voxelIter.pos()));
528 s -= target.wsSample(xyz);
529 s *= invMask ? 1 - a : a;
536 template<
typename Gr
idT,
typename InterruptT>
540 LevelSetMorphing<GridT, InterruptT>::
541 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
542 sampleAlignedSpeed(
const LeafRange& range,
Index speedBuffer)
544 typedef typename LeafType::ValueOnCIter VoxelIterT;
545 mParent->mTracker.checkInterrupter();
547 typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
549 if (mParent->mMask == NULL) {
550 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
551 BufferType& source = leafIter.buffer(speedBuffer);
552 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
553 ScalarType& s =
const_cast<ScalarType&
>(source.getValue(voxelIter.pos()));
554 s -= target.getValue(voxelIter.getCoord());
559 const ScalarType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
560 const bool invMask = mParent->isMaskInverted();
561 typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
562 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
563 BufferType& source = leafIter.buffer(speedBuffer);
564 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
565 const Coord ijk = voxelIter.getCoord();
567 ScalarType& s =
const_cast<ScalarType&
>(source.getValue(voxelIter.pos()));
568 s -= target.getValue(ijk);
569 s *= invMask ? 1 - a : a;
576 template<
typename Gr
idT,
typename InterruptT>
580 LevelSetMorphing<GridT, InterruptT>::
581 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
582 cook(ThreadingMode mode,
size_t swapBuffer)
584 mParent->mTracker.startInterrupter(
"Morphing level set");
586 const int grainSize = mParent->mTracker.getGrainSize();
587 const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
589 if (mParent->mTracker.getGrainSize()==0) {
591 }
else if (mode == PARALLEL_FOR) {
592 tbb::parallel_for(range, *
this);
593 }
else if (mode == PARALLEL_REDUCE) {
594 tbb::parallel_reduce(range, *
this);
596 throw std::runtime_error(
"Undefined threading mode");
599 mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
601 mParent->mTracker.endInterrupter();
606 template<
typename Gr
idT,
typename InterruptT>
607 template <
typename MapT,
611 LevelSetMorphing<GridT, InterruptT>::
612 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
613 euler1(
const LeafRange& range, ScalarType dt,
Index resultBuffer,
Index speedBuffer)
615 typedef math::BIAS_SCHEME<SpatialScheme> SchemeT;
616 typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
617 typedef typename LeafType::ValueOnCIter VoxelIterT;
618 typedef math::GradientNormSqrd<MapT, SpatialScheme> NumGrad;
620 mParent->mTracker.checkInterrupter();
621 const MapT& map = *mMap;
622 StencilT stencil(mParent->mTracker.grid());
624 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
625 BufferType& speed = leafIter.buffer(speedBuffer);
626 BufferType& result = leafIter.buffer(resultBuffer);
627 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
628 const Index n = voxelIter.pos();
629 const ScalarType S = speed.getValue(n);
631 stencil.moveTo(voxelIter);
632 result.setValue(n, *voxelIter - dt * S * NumGrad::result(map, stencil));
639 template<
typename Gr
idT,
typename InterruptT>
643 LevelSetMorphing<GridT, InterruptT>::
644 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
645 euler2(
const LeafRange& range, ScalarType dt, ScalarType alpha,
648 typedef math::BIAS_SCHEME<SpatialScheme> SchemeT;
649 typedef typename SchemeT::template ISStencil<GridType>::StencilType StencilT;
650 typedef typename LeafType::ValueOnCIter VoxelIterT;
651 typedef math::GradientNormSqrd<MapT, SpatialScheme> NumGrad;
653 mParent->mTracker.checkInterrupter();
654 const MapT& map = *mMap;
655 const ScalarType beta = ScalarType(1.0) - alpha;
656 StencilT stencil(mParent->mTracker.grid());
658 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
659 BufferType& speed = leafIter.buffer(speedBuffer);
660 BufferType& result = leafIter.buffer(resultBuffer);
661 BufferType& phi = leafIter.buffer(phiBuffer);
662 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
663 const Index n = voxelIter.pos();
664 const ScalarType S = speed.getValue(n);
666 stencil.moveTo(voxelIter);
667 const ScalarType G = NumGrad::result(map, stencil);
668 result.setValue(n, alpha * phi.getValue(n) + beta * (*voxelIter - dt * S * G));
677 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
Definition: LeafManager.h:126
Definition: FiniteDifference.h:265
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:284
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:606
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:693
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Index32 Index
Definition: Types.h:59
math::Vec3< Real > Vec3R
Definition: Types.h:77
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:263
Definition: Exceptions.h:39
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:326
Definition: FiniteDifference.h:264
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:545
Definition: Exceptions.h:88
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Definition: FiniteDifference.h:198
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Definition: FiniteDifference.h:194
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:314
Definition: FiniteDifference.h:263
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261