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>
86 LevelSetMorphing(GridT& sourceGrid,
const GridT& targetGrid, InterruptT* interrupt = NULL):
87 mTracker(sourceGrid, interrupt), mTarget(&targetGrid),
89 mTemporalScheme(math::
TVD_RK2) {}
94 void setTarget(
const GridT& targetGrid) { mTarget = &targetGrid; }
109 return mTracker.getSpatialScheme();
114 mTracker.setSpatialScheme(scheme);
119 return mTracker.getTemporalScheme();
124 mTracker.setTemporalScheme(scheme);
141 size_t advect(ScalarType time0, ScalarType time1);
152 LevelSetMorph(TrackerT& tracker,
const GridT* target);
154 LevelSetMorph(
const LevelSetMorph& other);
156 LevelSetMorph(LevelSetMorph& other, tbb::split);
158 virtual ~LevelSetMorph() {}
161 size_t advect(ScalarType time0, ScalarType time1);
163 void operator()(
const LeafRange& r)
const
165 if (mTask) mTask(const_cast<LevelSetMorph*>(
this), r);
169 void operator()(
const LeafRange& r)
171 if (mTask) mTask(
this, r);
175 void join(
const LevelSetMorph& other) { mMaxAbsS =
math::Max(mMaxAbsS, other.mMaxAbsS); }
177 typedef typename boost::function<void (LevelSetMorph*, const LeafRange&)> FuncType;
179 const GridT* mTarget;
180 ScalarType mMinAbsS, mMaxAbsS;
185 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
187 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
190 typename GridT::ValueType sampleSpeed(
191 ScalarType time0, ScalarType time1,
Index speedBuffer);
192 void sampleXformedSpeed(
const LeafRange& r,
Index speedBuffer);
193 void sampleAlignedSpeed(
const LeafRange& r,
Index speedBuffer);
196 void euler1(
const LeafRange& r, ScalarType dt,
Index resultBuffer,
Index speedBuffer);
200 void euler2(
const LeafRange& r, ScalarType dt, ScalarType alpha,
205 template<math::BiasedGradientScheme SpatialScheme>
206 size_t advect1(ScalarType time0, ScalarType time1);
210 size_t advect2(ScalarType time0, ScalarType time1);
215 size_t advect3(ScalarType time0, ScalarType time1);
218 const GridT* mTarget;
223 void operator=(
const LevelSetMorphing& other) {}
227 template<
typename Gr
idT,
typename InterruptT>
231 switch (mSpatialScheme) {
233 return this->advect1<math::FIRST_BIAS >(time0, time1);
241 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
248 template<
typename Gr
idT,
typename InterruptT>
249 template<math::BiasedGradientScheme SpatialScheme>
253 switch (mTemporalScheme) {
255 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
257 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
259 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
266 template<
typename Gr
idT,
typename InterruptT>
270 LevelSetMorphing<GridT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
273 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
274 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
275 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
276 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
278 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
279 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
280 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
281 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
288 template<
typename Gr
idT,
typename InterruptT>
293 LevelSetMorphing<GridT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
295 LevelSetMorph<MapT, SpatialScheme, TemporalScheme> tmp(mTracker, mTarget);
296 return tmp.advect(time0, time1);
303 template<
typename Gr
idT,
typename InterruptT>
307 LevelSetMorphing<GridT, InterruptT>::
308 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
309 LevelSetMorph(TrackerT& tracker,
const GridT* target):
313 mMap(tracker.grid().transform().template constMap<MapT>().get()),
318 template<
typename Gr
idT,
typename InterruptT>
322 LevelSetMorphing<GridT, InterruptT>::
323 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
324 LevelSetMorph(
const LevelSetMorph& other):
325 mTracker(other.mTracker),
326 mTarget(other.mTarget),
327 mMinAbsS(other.mMinAbsS),
328 mMaxAbsS(other.mMaxAbsS),
334 template<
typename Gr
idT,
typename InterruptT>
338 LevelSetMorphing<GridT, InterruptT>::
339 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
340 LevelSetMorph(LevelSetMorph& other, tbb::split):
341 mTracker(other.mTracker),
342 mTarget(other.mTarget),
343 mMinAbsS(other.mMinAbsS),
344 mMaxAbsS(other.mMaxAbsS),
350 template<
typename Gr
idT,
typename InterruptT>
356 advect(ScalarType time0, ScalarType time1)
362 while (time0 < time1 && mTracker->checkInterrupter()) {
363 mTracker->leafs().rebuildAuxBuffers(auxBuffers);
365 const ScalarType dt = this->sampleSpeed(time0, time1, auxBuffers);
369 switch(TemporalScheme) {
373 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 2);
375 this->cook(PARALLEL_FOR, 1);
380 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 2);
382 this->cook(PARALLEL_FOR, 1);
386 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.5),
389 this->cook(PARALLEL_FOR, 1);
394 mTask = boost::bind(&LevelSetMorph::euler1, _1, _2, dt, 1, 3);
396 this->cook(PARALLEL_FOR, 1);
400 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(0.75),
403 this->cook(PARALLEL_FOR, 2);
407 mTask = boost::bind(&LevelSetMorph::euler2, _1, _2, dt, ScalarType(1.0/3.0),
410 this->cook(PARALLEL_FOR, 2);
413 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
419 mTracker->leafs().removeAuxBuffers();
428 template<
typename Gr
idT,
typename InterruptT>
431 inline typename GridT::ValueType
432 LevelSetMorphing<GridT, InterruptT>::
433 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
434 sampleSpeed(ScalarType time0, ScalarType time1,
Index speedBuffer)
437 const size_t leafCount = mTracker->leafs().leafCount();
438 if (leafCount==0 || time0 >= time1)
return ScalarType(0);
440 if (mTarget->transform() == mTracker->grid().transform()) {
441 mTask = boost::bind(&LevelSetMorph::sampleAlignedSpeed, _1, _2, speedBuffer);
443 mTask = boost::bind(&LevelSetMorph::sampleXformedSpeed, _1, _2, speedBuffer);
445 this->cook(PARALLEL_REDUCE);
447 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
448 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
449 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
450 const ScalarType dt =
math::Abs(time1 - time0), dx = mTracker->voxelSize();
451 return math::Min(dt, ScalarType(CFL*dx/mMaxAbsS));
454 template<
typename Gr
idT,
typename InterruptT>
458 LevelSetMorphing<GridT, InterruptT>::
459 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
460 sampleXformedSpeed(
const LeafRange& range,
Index speedBuffer)
462 typedef typename LeafType::ValueOnCIter VoxelIterT;
463 typedef tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler> SamplerT;
464 const MapT& map = *mMap;
465 mTracker->checkInterrupter();
467 SamplerT sampler(mTarget->getAccessor(), mTarget->transform());
468 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
469 BufferType& speed = leafIter.buffer(speedBuffer);
470 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
471 ScalarType& s =
const_cast<ScalarType&
>(speed.getValue(voxelIter.pos()));
472 s -= sampler.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
478 template<
typename Gr
idT,
typename InterruptT>
482 LevelSetMorphing<GridT, InterruptT>::
483 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
484 sampleAlignedSpeed(
const LeafRange& range,
Index speedBuffer)
486 typedef typename LeafType::ValueOnCIter VoxelIterT;
487 mTracker->checkInterrupter();
489 typename GridT::ConstAccessor target = mTarget->getAccessor();
490 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
491 BufferType& speed = leafIter.buffer(speedBuffer);
492 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
493 ScalarType& s =
const_cast<ScalarType&
>(speed.getValue(voxelIter.pos()));
494 s -= target.getValue(voxelIter.getCoord());
500 template<
typename Gr
idT,
typename InterruptT>
504 LevelSetMorphing<GridT, InterruptT>::
505 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
506 cook(ThreadingMode mode,
size_t swapBuffer)
508 mTracker->startInterrupter(
"Morphing level set");
510 const int grainSize = mTracker->getGrainSize();
511 const LeafRange range = mTracker->leafs().leafRange(grainSize);
513 if (mTracker->getGrainSize()==0) {
515 }
else if (mode == PARALLEL_FOR) {
516 tbb::parallel_for(range, *
this);
517 }
else if (mode == PARALLEL_REDUCE) {
518 tbb::parallel_reduce(range, *
this);
520 throw std::runtime_error(
"Undefined threading mode");
523 mTracker->leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
525 mTracker->endInterrupter();
530 template<
typename Gr
idT,
typename InterruptT>
531 template <
typename MapT,
535 LevelSetMorphing<GridT, InterruptT>::
536 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
537 euler1(
const LeafRange& range, ScalarType dt,
Index resultBuffer,
Index speedBuffer)
539 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
540 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
541 typedef typename LeafType::ValueOnCIter VoxelIterT;
543 mTracker->checkInterrupter();
544 const MapT& map = *mMap;
545 Stencil stencil(mTracker->grid());
547 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
548 BufferType& speed = leafIter.buffer(speedBuffer);
549 BufferType& result = leafIter.buffer(resultBuffer);
550 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
551 const Index n = voxelIter.pos();
552 stencil.moveTo(voxelIter);
554 result.setValue(n, *voxelIter - dt * speed.getValue(n) * G);
561 template<
typename Gr
idT,
typename InterruptT>
565 LevelSetMorphing<GridT, InterruptT>::
566 LevelSetMorph<MapT, SpatialScheme, TemporalScheme>::
567 euler2(
const LeafRange& range, ScalarType dt, ScalarType alpha,
570 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
571 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
572 typedef typename LeafType::ValueOnCIter VoxelIterT;
574 mTracker->checkInterrupter();
575 const MapT& map = *mMap;
576 const ScalarType beta = ScalarType(1.0) - alpha;
577 Stencil stencil(mTracker->grid());
579 for (
typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
580 BufferType& speed = leafIter.buffer(speedBuffer);
581 BufferType& result = leafIter.buffer(resultBuffer);
582 BufferType& phi = leafIter.buffer(phiBuffer);
583 for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
584 const Index n = voxelIter.pos();
585 stencil.moveTo(voxelIter);
588 alpha * phi.getValue(n) + beta * (*voxelIter - dt * speed.getValue(n) * G));
597 #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
static Accessor::ValueType result(const MapType &map, const Accessor &grid, const Coord &ijk)
Definition: Operators.h:848
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:653
Index32 Index
Definition: Types.h:57
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:244
Definition: FiniteDifference.h:264
Definition: Exceptions.h:88
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:505
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
Definition: LeafManager.h:122
Definition: FiniteDifference.h:265
#define OPENVDB_VERSION_NAME
Definition: version.h:45
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Definition: FiniteDifference.h:198
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
Definition: FiniteDifference.h:263
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:274
Definition: FiniteDifference.h:194
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:566