11 #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
12 #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
45 template<
typename GridT,
46 typename InterruptT = util::NullInterrupter>
59 LevelSetMorphing(GridT& sourceGrid,
const GridT& targetGrid, InterruptT* interrupt =
nullptr)
60 : mTracker(sourceGrid, interrupt)
61 , mTarget(&targetGrid)
64 , mTemporalScheme(math::
TVD_RK2)
74 void setTarget(
const GridT& targetGrid) { mTarget = &targetGrid; }
137 if (!(min < max))
OPENVDB_THROW(ValueError,
"Invalid mask range (expects min < max)");
139 mDeltaMask = max-
min;
161 template<math::BiasedGradientScheme SpatialScheme>
174 const GridT *mTarget, *mMask;
188 Morph(
const Morph& other);
199 if (mTask) mTask(const_cast<Morph*>(
this), r);
200 else OPENVDB_THROW(ValueError,
"task is undefined - don\'t call this method directly");
205 if (mTask) mTask(
this, r);
206 else OPENVDB_THROW(ValueError,
"task is undefined - don\'t call this method directly");
209 void join(
const Morph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
212 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
214 void cook(ThreadingMode
mode,
size_t swapBuffer = 0);
218 void sampleXformedSpeed(
const LeafRange& r, Index speedBuffer);
219 void sampleAlignedSpeed(
const LeafRange& r, Index speedBuffer);
223 template <
int Nominator,
int Denominator>
230 using FuncType =
typename std::function<void (Morph*, const LeafRange&)>;
239 template<
typename Gr
idT,
typename InterruptT>
243 switch (mSpatialScheme) {
245 return this->advect1<math::FIRST_BIAS >(time0, time1);
253 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
259 OPENVDB_THROW(ValueError,
"Spatial difference scheme not supported!");
264 template<
typename Gr
idT,
typename InterruptT>
265 template<math::BiasedGradientScheme SpatialScheme>
269 switch (mTemporalScheme) {
271 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
273 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
275 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
278 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
283 template<
typename Gr
idT,
typename InterruptT>
289 const math::Transform&
trans = mTracker.grid().transform();
291 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
293 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
296 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
298 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
305 template<
typename Gr
idT,
typename InterruptT>
312 Morph<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
313 return tmp.advect(time0, time1);
319 template<
typename Gr
idT,
typename InterruptT>
323 LevelSetMorphing<GridT, InterruptT>::
324 Morph<MapT, SpatialScheme, TemporalScheme>::
325 Morph(LevelSetMorphing<GridT, InterruptT>& parent)
328 , mMap(parent.mTracker.grid().
transform().template constMap<MapT>().
get())
333 template<
typename Gr
idT,
typename InterruptT>
337 LevelSetMorphing<GridT, InterruptT>::
338 Morph<MapT, SpatialScheme, TemporalScheme>::
339 Morph(
const Morph& other)
340 : mParent(other.mParent)
341 , mMinAbsS(other.mMinAbsS)
342 , mMaxAbsS(other.mMaxAbsS)
348 template<
typename Gr
idT,
typename InterruptT>
352 LevelSetMorphing<GridT, InterruptT>::
353 Morph<MapT, SpatialScheme, TemporalScheme>::
355 : mParent(other.mParent)
356 , mMinAbsS(other.mMinAbsS)
357 , mMaxAbsS(other.mMaxAbsS)
363 template<
typename Gr
idT,
typename InterruptT>
371 namespace ph = std::placeholders;
377 while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
378 mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
380 const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
384 switch(TemporalScheme) {
388 mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, 2);
391 this->cook(PARALLEL_FOR, 1);
396 mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, 2);
399 this->cook(PARALLEL_FOR, 1);
403 mTask = std::bind(&Morph::euler12, ph::_1, ph::_2, dt);
406 this->cook(PARALLEL_FOR, 1);
411 mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, 3);
414 this->cook(PARALLEL_FOR, 1);
418 mTask = std::bind(&Morph::euler34, ph::_1, ph::_2, dt);
421 this->cook(PARALLEL_FOR, 2);
425 mTask = std::bind(&Morph::euler13, ph::_1, ph::_2, dt);
428 this->cook(PARALLEL_FOR, 2);
432 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
438 mParent->mTracker.leafs().removeAuxBuffers();
441 mParent->mTracker.track();
447 template<
typename Gr
idT,
typename InterruptT>
451 LevelSetMorphing<GridT, InterruptT>::
452 Morph<MapT, SpatialScheme, TemporalScheme>::
455 namespace ph = std::placeholders;
458 const size_t leafCount = mParent->mTracker.leafs().leafCount();
459 if (leafCount==0 || time0 >= time1)
return ValueType(0);
461 const math::Transform& xform = mParent->mTracker.grid().transform();
462 if (mParent->mTarget->transform() == xform &&
463 (mParent->mMask ==
nullptr || mParent->mMask->transform() == xform)) {
464 mTask = std::bind(&Morph::sampleAlignedSpeed, ph::_1, ph::_2, speedBuffer);
466 mTask = std::bind(&Morph::sampleXformedSpeed, ph::_1, ph::_2, speedBuffer);
468 this->cook(PARALLEL_REDUCE);
471 TemporalScheme == math::
TVD_RK2 ? ValueType(0.9) :
472 ValueType(1.0))/math::
Sqrt(ValueType(3.0));
473 const ValueType dt =
math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
477 template<
typename Gr
idT,
typename InterruptT>
481 LevelSetMorphing<GridT, InterruptT>::
482 Morph<MapT, SpatialScheme, TemporalScheme>::
483 sampleXformedSpeed(
const LeafRange&
range,
Index speedBuffer)
485 using VoxelIterT =
typename LeafType::ValueOnCIter;
486 using SamplerT = tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler>;
488 const MapT& map = *mMap;
489 mParent->mTracker.checkInterrupter();
491 typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
492 SamplerT
target(targetAcc, mParent->mTarget->transform());
493 if (mParent->mMask ==
nullptr) {
494 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
495 ValueType* speed = leafIter.buffer(speedBuffer).data();
497 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
498 ValueType&
s = speed[voxelIter.pos()];
499 s -=
target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
506 const ValueType
min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
507 const bool invMask = mParent->isMaskInverted();
508 typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
509 SamplerT
mask(maskAcc, mParent->mMask->transform());
510 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
511 ValueType* speed = leafIter.buffer(speedBuffer).data();
513 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
514 const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());
516 ValueType& s = speed[voxelIter.pos()];
517 s -=
target.wsSample(xyz);
518 s *= invMask ? 1 - a :
a;
527 template<
typename Gr
idT,
typename InterruptT>
531 LevelSetMorphing<GridT, InterruptT>::
532 Morph<MapT, SpatialScheme, TemporalScheme>::
533 sampleAlignedSpeed(
const LeafRange& range,
Index speedBuffer)
535 using VoxelIterT =
typename LeafType::ValueOnCIter;
537 mParent->mTracker.checkInterrupter();
539 typename GridT::ConstAccessor
target = mParent->mTarget->getAccessor();
541 if (mParent->mMask ==
nullptr) {
542 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
543 ValueType* speed = leafIter.buffer(speedBuffer).data();
545 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
546 ValueType& s = speed[voxelIter.pos()];
547 s -= target.getValue(voxelIter.getCoord());
554 const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
555 const bool invMask = mParent->isMaskInverted();
556 typename GridT::ConstAccessor
mask = mParent->mMask->getAccessor();
557 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
558 ValueType* speed = leafIter.buffer(speedBuffer).data();
560 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
561 const Coord ijk = voxelIter.getCoord();
563 ValueType& s = speed[voxelIter.pos()];
564 s -= target.getValue(ijk);
565 s *= invMask ? 1 - a :
a;
574 template<
typename Gr
idT,
typename InterruptT>
578 LevelSetMorphing<GridT, InterruptT>::
579 Morph<MapT, SpatialScheme, TemporalScheme>::
580 cook(ThreadingMode
mode,
size_t swapBuffer)
582 mParent->mTracker.startInterrupter(
"Morphing level set");
584 const int grainSize = mParent->mTracker.getGrainSize();
585 const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
587 if (mParent->mTracker.getGrainSize()==0) {
589 }
else if (mode == PARALLEL_FOR) {
591 }
else if (mode == PARALLEL_REDUCE) {
592 tbb::parallel_reduce(range, *
this);
594 OPENVDB_THROW(ValueError,
"expected threading mode " <<
int(PARALLEL_FOR)
595 <<
" or " <<
int(PARALLEL_REDUCE) <<
", got " <<
int(mode));
598 mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
600 mParent->mTracker.endInterrupter();
603 template<
typename Gr
idT,
typename InterruptT>
606 template <
int Nominator,
int Denominator>
608 LevelSetMorphing<GridT,InterruptT>::
609 Morph<MapT, SpatialScheme, TemporalScheme>::
610 euler(
const LeafRange& range, ValueType dt,
613 using SchemeT = math::BIAS_SCHEME<SpatialScheme>;
614 using StencilT =
typename SchemeT::template ISStencil<GridType>::StencilType;
615 using VoxelIterT =
typename LeafType::ValueOnCIter;
616 using NumGrad = math::GradientNormSqrd<MapT, SpatialScheme>;
621 mParent->mTracker.checkInterrupter();
622 const MapT& map = *mMap;
623 StencilT
stencil(mParent->mTracker.grid());
625 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
626 const ValueType* speed = leafIter.buffer(speedBuffer).data();
628 const ValueType* phi = leafIter.buffer(phiBuffer).data();
629 ValueType*
result = leafIter.buffer(resultBuffer).data();
630 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
631 const Index n = voxelIter.pos();
635 result[
n] = Nominator ? Alpha * phi[
n] + Beta * v :
v;
646 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
648 #ifdef OPENVDB_INSTANTIATE_LEVELSETMORPH
655 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
662 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
TemporalIntegrationScheme
Temporal integration schemes.
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
GLboolean GLboolean GLboolean GLboolean a
#define OPENVDB_USE_VERSION_NAMESPACE
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
**But if you need a result
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
#define OPENVDB_INSTANTIATE_CLASS
Coord Abs(const Coord &xyz)
GA_API const UT_StringHolder trans
auto get(const UT_ARTIterator< T > &it) -> decltype(it.key())
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
GA_API const UT_StringHolder transform
static Name mapType()
Return UnitaryMap.
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
GLint GLfloat GLint stencil
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x .
__hostdev__ float Sqrt(float x)
Return the square root of a floating-point value.
GA_API const UT_StringHolder Alpha
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
#define OPENVDB_THROW(exception, message)