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>
80 template<
typename GridT,
81 typename InterruptT = util::NullInterrupter>
94 LevelSetMorphing(GridT& sourceGrid,
const GridT& targetGrid, InterruptT* interrupt = NULL):
95 mTracker(sourceGrid, interrupt), mTarget(&targetGrid),
97 mTemporalScheme(math::
TVD_RK2) {}
102 void setTarget(
const GridT& targetGrid) { mTarget = &targetGrid; }
117 return mTracker.getSpatialScheme();
122 mTracker.setSpatialScheme(scheme);
127 return mTracker.getTemporalScheme();
132 mTracker.setTemporalScheme(scheme);
149 size_t advect(ScalarType time0, ScalarType time1);
160 LevelSetMorph(TrackerT& tracker,
const GridT* target);
162 LevelSetMorph(
const LevelSetMorph& other);
164 LevelSetMorph(LevelSetMorph& other, tbb::split);
166 virtual ~LevelSetMorph() {}
169 size_t advect(ScalarType time0, ScalarType time1);
171 void operator()(
const LeafRange& r)
const
173 if (mTask) mTask(const_cast<LevelSetMorph*>(
this), r);
177 void operator()(
const LeafRange& r)
179 if (mTask) mTask(
this, r);
183 void join(
const LevelSetMorph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
185 typedef typename boost::function<void (LevelSetMorph*, const LeafRange&)> FuncType;
187 const GridT* mTarget;
188 ScalarType mMinAbsS, mMaxAbsS;
193 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
195 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
198 typename GridT::ValueType sampleSpeed(
199 ScalarType time0, ScalarType time1,
Index speedBuffer);
200 void sampleXformedSpeed(
const LeafRange& r,
Index speedBuffer);
201 void sampleAlignedSpeed(
const LeafRange& r,
Index speedBuffer);
204 void euler1(
const LeafRange& r, ScalarType dt,
Index resultBuffer,
Index speedBuffer);
208 void euler2(
const LeafRange& r, ScalarType dt, ScalarType alpha,
213 template<math::BiasedGradientScheme SpatialScheme>
214 size_t advect1(ScalarType time0, ScalarType time1);
218 size_t advect2(ScalarType time0, ScalarType time1);
223 size_t advect3(ScalarType time0, ScalarType time1);
226 const GridT* mTarget;
231 void operator=(
const LevelSetMorphing& other) {}
235 template<
typename Gr
idT,
typename InterruptT>
239 switch (mSpatialScheme) {
241 return this->advect1<math::FIRST_BIAS >(time0, time1);
243 return this->advect1<math::SECOND_BIAS >(time0, time1);
245 return this->advect1<math::THIRD_BIAS >(time0, time1);
247 return this->advect1<math::WENO5_BIAS >(time0, time1);
249 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
256 template<
typename Gr
idT,
typename InterruptT>
257 template<math::BiasedGradientScheme SpatialScheme>
261 switch (mTemporalScheme) {
263 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
265 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
267 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
274 template<
typename Gr
idT,
typename InterruptT>
278 LevelSetMorphing<GridT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
281 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
282 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
283 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
284 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
286 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
287 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
288 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
289 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
296 template<
typename Gr
idT,
typename InterruptT>
301 LevelSetMorphing<GridT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
303 LevelSetMorph<MapT, SpatialScheme, TemporalScheme> tmp(mTracker, mTarget);
304 return tmp.advect(time0, time1);
311 template<
typename Gr
idT,
typename InterruptT>
315 LevelSetMorphing<GridT, InterruptT>::
316 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
317 LevelSetMorph(TrackerT& tracker,
const GridT* target):
321 mMap(tracker.grid().transform().template constMap<MapT>().get()),
326 template<
typename Gr
idT,
typename InterruptT>
330 LevelSetMorphing<GridT, InterruptT>::
331 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
332 LevelSetMorph(
const LevelSetMorph& other):
333 mTracker(other.mTracker),
334 mTarget(other.mTarget),
335 mMinAbsS(other.mMinAbsS),
336 mMaxAbsS(other.mMaxAbsS),
342 template<
typename Gr
idT,
typename InterruptT>
346 LevelSetMorphing<GridT, InterruptT>::
347 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
348 LevelSetMorph(LevelSetMorph& other, tbb::split):
349 mTracker(other.mTracker),
350 mTarget(other.mTarget),
351 mMinAbsS(other.mMinAbsS),
352 mMaxAbsS(other.mMaxAbsS),
358 template<
typename Gr
idT,
typename InterruptT>
364 advect(ScalarType time0, ScalarType time1)
370 while (time0 < time1 && mTracker->checkInterrupter()) {
371 mTracker->leafs().rebuildAuxBuffers(auxBuffers);
373 const ScalarType dt = this->sampleSpeed(time0, time1, auxBuffers);
377 switch(TemporalScheme) {
381 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 2);
383 this->cook(PARALLEL_FOR, 1);
388 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 2);
390 this->cook(PARALLEL_FOR, 1);
394 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.5),
397 this->cook(PARALLEL_FOR, 1);
402 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 3);
404 this->cook(PARALLEL_FOR, 1);
408 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.75),
411 this->cook(PARALLEL_FOR, 2);
415 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(1.0/3.0),
418 this->cook(PARALLEL_FOR, 2);
421 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
427 mTracker->leafs().removeAuxBuffers();
436 template<
typename Gr
idT,
typename InterruptT>
439 inline typename GridT::ValueType
440 LevelSetMorphing<GridT, InterruptT>::
441 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
442 sampleSpeed(ScalarType time0, ScalarType time1,
Index speedBuffer)
445 const size_t leafCount = mTracker->leafs().leafCount();
446 if (leafCount==0 || time0 >= time1)
return ScalarType(0);
448 if (mTarget->transform() == mTracker->grid().transform()) {
449 mTask = boost::bind(&LevelSetMorph::sampleAlignedSpeed, _1, _2, speedBuffer);
451 mTask = boost::bind(&LevelSetMorph::sampleXformedSpeed, _1, _2, speedBuffer);
453 this->cook(PARALLEL_REDUCE);
455 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
456 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
457 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
458 const ScalarType dt =
math::Abs(time1 - time0), dx = mTracker->voxelSize();
459 return math::Min(dt, ScalarType(CFL*dx/mMaxAbsS));
462 template<
typename Gr
idT,
typename InterruptT>
466 LevelSetMorphing<GridT, InterruptT>::
467 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
468 sampleXformedSpeed(
const LeafRange& range,
Index speedBuffer)
470 typedef typename LeafType::ValueOnCIter VoxelIterT;
471 typedef tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler> SamplerT;
472 const MapT& map = *
mMap;
473 mTracker->checkInterrupter();
475 SamplerT sampler(mTarget->getAccessor(), mTarget->transform());
476 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
477 BufferType& speed = leafIter.buffer(speedBuffer);
478 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
479 ScalarType& s =
const_cast<ScalarType&
>(speed.getValue(voxelIter.pos()));
480 s -= sampler.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
486 template<
typename Gr
idT,
typename InterruptT>
490 LevelSetMorphing<GridT, InterruptT>::
491 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
492 sampleAlignedSpeed(
const LeafRange& range,
Index speedBuffer)
494 typedef typename LeafType::ValueOnCIter VoxelIterT;
495 mTracker->checkInterrupter();
497 typename GridT::ConstAccessor target = mTarget->getAccessor();
498 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
499 BufferType& speed = leafIter.buffer(speedBuffer);
500 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
501 ScalarType& s =
const_cast<ScalarType&
>(speed.getValue(voxelIter.pos()));
502 s -= target.getValue(voxelIter.getCoord());
508 template<
typename Gr
idT,
typename InterruptT>
512 LevelSetMorphing<GridT, InterruptT>::
513 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
514 cook(ThreadingMode mode,
size_t swapBuffer)
516 mTracker->startInterrupter(
"Morphing level set");
518 const int grainSize = mTracker->getGrainSize();
519 const LeafRange range = mTracker->leafs().leafRange(grainSize);
521 if (mTracker->getGrainSize()==0) {
523 }
else if (mode == PARALLEL_FOR) {
524 tbb::parallel_for(range, *
this);
525 }
else if (mode == PARALLEL_REDUCE) {
526 tbb::parallel_reduce(range, *
this);
528 throw std::runtime_error(
"Undefined threading mode");
531 mTracker->leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
533 mTracker->endInterrupter();
538 template<
typename Gr
idT,
typename InterruptT>
539 template <
typename MapT,
543 LevelSetMorphing<GridT, InterruptT>::
544 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
545 euler1(
const LeafRange& range, ScalarType dt,
Index resultBuffer,
Index speedBuffer)
547 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
548 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
549 typedef typename LeafType::ValueOnCIter VoxelIterT;
551 mTracker->checkInterrupter();
552 const MapT& map = *
mMap;
553 Stencil stencil(mTracker->grid());
555 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
556 BufferType& speed = leafIter.buffer(speedBuffer);
557 BufferType& result = leafIter.buffer(resultBuffer);
558 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
559 const Index n = voxelIter.pos();
560 stencil.moveTo(voxelIter);
562 result.setValue(n, *voxelIter - dt * speed.getValue(n) * G);
569 template<
typename Gr
idT,
typename InterruptT>
573 LevelSetMorphing<GridT, InterruptT>::
574 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
575 euler2(
const LeafRange& range, ScalarType dt, ScalarType alpha,
578 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
579 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
580 typedef typename LeafType::ValueOnCIter VoxelIterT;
582 mTracker->checkInterrupter();
583 const MapT& map = *
mMap;
584 const ScalarType beta = ScalarType(1.0) - alpha;
585 Stencil stencil(mTracker->grid());
587 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
588 BufferType& speed = leafIter.buffer(speedBuffer);
589 BufferType& result = leafIter.buffer(resultBuffer);
590 BufferType& phi = leafIter.buffer(phiBuffer);
591 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
592 const Index n = voxelIter.pos();
593 stencil.moveTo(voxelIter);
596 alpha * phi.getValue(n) + beta * (*voxelIter - dt * speed.getValue(n) * G));
605 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
Definition: FiniteDifference.h:196
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
TemporalIntegrationScheme
Temporal integrations schemes.
Definition: FiniteDifference.h:261
const MapT & mMap
Definition: GridOperators.h:384
Definition: FiniteDifference.h:265
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Definition: Exceptions.h:88
Definition: FiniteDifference.h:264
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:253
Definition: FiniteDifference.h:263
Index32 Index
Definition: Types.h:56
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:655
Definition: FiniteDifference.h:198
#define OPENVDB_VERSION_NAME
Definition: version.h:45
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:575
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:283
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:836
Definition: LeafManager.h:122
Definition: FiniteDifference.h:195
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
Definition: FiniteDifference.h:197
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:514
Definition: FiniteDifference.h:194