Loading...
Searching...
No Matches
LevelSetMorph.h
Go to the documentation of this file.
// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0
/// @author Ken Museth
///
/// @file tools/LevelSetMorph.h
///
/// @brief Shape morphology of level sets. Morphing from a source
/// narrow-band level sets to a target narrow-band level set.
#ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
#define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
#include "LevelSetTracker.h"
#include "Interpolation.h" // for BoxSampler, etc.
#include <openvdb/math/FiniteDifference.h>
#include <functional>
#include <limits>
namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace tools {
/// @brief Shape morphology of level sets. Morphing from a source
/// narrow-band level sets to a target narrow-band level set.
///
/// @details
/// The @c InterruptType template argument below refers to any class
/// with the following interface:
/// @code
/// class Interrupter {
/// ...
/// public:
/// void start(const char* name = nullptr) // called when computations begin
/// void end() // called when computations end
/// bool wasInterrupted(int percent=-1) // return true to break computation
/// };
/// @endcode
///
/// @note If no template argument is provided for this InterruptType,
/// the util::NullInterrupter is used, which implies that all interrupter
/// calls are no-ops (i.e., they incur no computational overhead).
template<typename GridT,
typename InterruptT = util::NullInterrupter>
class LevelSetMorphing
{
public:
using GridType = GridT;
using TreeType = typename GridT::TreeType;
using TrackerT = LevelSetTracker<GridT, InterruptT>;
using LeafRange = typename TrackerT::LeafRange;
using LeafType = typename TrackerT::LeafType;
using BufferType = typename TrackerT::BufferType;
using ValueType = typename TrackerT::ValueType;
/// Main constructor
LevelSetMorphing(GridT& sourceGrid, const GridT& targetGrid, InterruptT* interrupt = nullptr)
: mTracker(sourceGrid, interrupt)
, mTarget(&targetGrid)
, mMask(nullptr)
, mSpatialScheme(math::HJWENO5_BIAS)
, mTemporalScheme(math::TVD_RK2)
, mMinMask(0)
, mDeltaMask(1)
, mInvertMask(false)
{
}
virtual ~LevelSetMorphing() {}
/// Redefine the target level set
void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; }
/// Define the alpha mask
void setAlphaMask(const GridT& maskGrid) { mMask = &maskGrid; }
/// Return the spatial finite-difference scheme
math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; }
/// Set the spatial finite-difference scheme
void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; }
/// Return the temporal integration scheme
math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; }
/// Set the temporal integration scheme
void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; }
/// Return the spatial finite-difference scheme
math::BiasedGradientScheme getTrackerSpatialScheme() const
{
return mTracker.getSpatialScheme();
}
/// Set the spatial finite-difference scheme
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
{
mTracker.setSpatialScheme(scheme);
}
/// Return the temporal integration scheme
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
{
return mTracker.getTemporalScheme();
}
/// Set the temporal integration scheme
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
{
mTracker.setTemporalScheme(scheme);
}
/// Return the number of normalizations performed per track or normalize call.
int getNormCount() const { return mTracker.getNormCount(); }
/// Set the number of normalizations performed per track or normalize call.
void setNormCount(int n) { mTracker.setNormCount(n); }
/// Return the grain size used for multithreading
int getGrainSize() const { return mTracker.getGrainSize(); }
/// @brief Set the grain size used for multithreading.
/// @note A grain size of 0 or less disables multithreading!
void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); }
/// @brief Return the minimum value of the mask to be used for the
/// derivation of a smooth alpha value.
ValueType minMask() const { return mMinMask; }
/// @brief Return the maximum value of the mask to be used for the
/// derivation of a smooth alpha value.
ValueType maxMask() const { return mDeltaMask + mMinMask; }
/// @brief Define the range for the (optional) scalar mask.
/// @param min Minimum value of the range.
/// @param max Maximum value of the range.
/// @details Mask values outside the range maps to alpha values of
/// respectfully zero and one, and values inside the range maps
/// smoothly to 0->1 (unless of course the mask is inverted).
/// @throw ValueError if @a min is not smaller than @a max.
void setMaskRange(ValueType min, ValueType max)
{
if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
mMinMask = min;
mDeltaMask = max-min;
}
/// @brief Return true if the mask is inverted, i.e. min->max in the
/// original mask maps to 1->0 in the inverted alpha mask.
bool isMaskInverted() const { return mInvertMask; }
/// @brief Invert the optional mask, i.e. min->max in the original
/// mask maps to 1->0 in the inverted alpha mask.
void invertMask(bool invert=true) { mInvertMask = invert; }
/// @brief Advect the level set from its current time, @a time0, to its
/// final time, @a time1. If @a time0 > @a time1, perform backward advection.
///
/// @return the number of CFL iterations used to advect from @a time0 to @a time1
size_t advect(ValueType time0, ValueType time1);
private:
// disallow copy construction and copy by assignment!
LevelSetMorphing(const LevelSetMorphing&);// not implemented
LevelSetMorphing& operator=(const LevelSetMorphing&);// not implemented
template<math::BiasedGradientScheme SpatialScheme>
size_t advect1(ValueType time0, ValueType time1);
template<math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
size_t advect2(ValueType time0, ValueType time1);
template<math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme,
typename MapType>
size_t advect3(ValueType time0, ValueType time1);
TrackerT mTracker;
const GridT *mTarget, *mMask;
math::BiasedGradientScheme mSpatialScheme;
math::TemporalIntegrationScheme mTemporalScheme;
ValueType mMinMask, mDeltaMask;
bool mInvertMask;
// This templated private class implements all the level set magic.
template<typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
struct Morph
{
/// Main constructor
Morph(LevelSetMorphing<GridT, InterruptT>& parent);
/// Shallow copy constructor called by tbb::parallel_for() threads
Morph(const Morph& other);
/// Shallow copy constructor called by tbb::parallel_reduce() threads
Morph(Morph& other, tbb::split);
/// destructor
virtual ~Morph() {}
/// Advect the level set from its current time, time0, to its final time, time1.
/// @return number of CFL iterations
size_t advect(ValueType time0, ValueType time1);
/// Used internally by tbb::parallel_for()
void operator()(const LeafRange& r) const
{
if (mTask) mTask(const_cast<Morph*>(this), r);
else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
}
/// Used internally by tbb::parallel_reduce()
void operator()(const LeafRange& r)
{
if (mTask) mTask(this, r);
else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
}
/// This is only called by tbb::parallel_reduce() threads
void join(const Morph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); }
/// Enum to define the type of multithreading
enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use
// method calling tbb
void cook(ThreadingMode mode, size_t swapBuffer = 0);
/// Sample field and return the CFT time step
typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer);
void sampleXformedSpeed(const LeafRange& r, Index speedBuffer);
void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer);
// Convex combination of Phi and a forward Euler advection steps:
// Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|);
template <int Nominator, int Denominator>
void euler(const LeafRange&, ValueType, Index, Index, Index);
inline void euler01(const LeafRange& r, ValueType t, Index s) {this->euler<0,1>(r,t,0,1,s);}
inline void euler12(const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);}
inline void euler34(const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);}
inline void euler13(const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);}
using FuncType = typename std::function<void (Morph*, const LeafRange&)>;
LevelSetMorphing* mParent;
ValueType mMinAbsS, mMaxAbsS;
const MapT* mMap;
FuncType mTask;
}; // end of private Morph struct
};//end of LevelSetMorphing
template<typename GridT, typename InterruptT>
inline size_t
LevelSetMorphing<GridT, InterruptT>::advect(ValueType time0, ValueType time1)
{
switch (mSpatialScheme) {
case math::FIRST_BIAS:
return this->advect1<math::FIRST_BIAS >(time0, time1);
//case math::SECOND_BIAS:
//return this->advect1<math::SECOND_BIAS >(time0, time1);
//case math::THIRD_BIAS:
//return this->advect1<math::THIRD_BIAS >(time0, time1);
//case math::WENO5_BIAS:
//return this->advect1<math::WENO5_BIAS >(time0, time1);
case math::HJWENO5_BIAS:
return this->advect1<math::HJWENO5_BIAS>(time0, time1);
case math::SECOND_BIAS:
case math::THIRD_BIAS:
case math::WENO5_BIAS:
case math::UNKNOWN_BIAS:
default:
OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
}
return 0;
}
template<typename GridT, typename InterruptT>
template<math::BiasedGradientScheme SpatialScheme>
inline size_t
LevelSetMorphing<GridT, InterruptT>::advect1(ValueType time0, ValueType time1)
{
switch (mTemporalScheme) {
case math::TVD_RK1:
return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
case math::TVD_RK2:
return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
case math::TVD_RK3:
return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
case math::UNKNOWN_TIS:
default:
OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
}
return 0;
}
template<typename GridT, typename InterruptT>
template<math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline size_t
LevelSetMorphing<GridT, InterruptT>::advect2(ValueType time0, ValueType time1)
{
const math::Transform& trans = mTracker.grid().transform();
if (trans.mapType() == math::UniformScaleMap::mapType()) {
return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
} else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) {
return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(
time0, time1);
} else if (trans.mapType() == math::UnitaryMap::mapType()) {
return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
} else if (trans.mapType() == math::TranslationMap::mapType()) {
return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
} else {
OPENVDB_THROW(ValueError, "MapType not supported!");
}
return 0;
}
template<typename GridT, typename InterruptT>
template<math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme,
typename MapT>
inline size_t
LevelSetMorphing<GridT, InterruptT>::advect3(ValueType time0, ValueType time1)
{
Morph<MapT, SpatialScheme, TemporalScheme> tmp(*this);
return tmp.advect(time0, time1);
}
///////////////////////////////////////////////////////////////////////
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
Morph(LevelSetMorphing<GridT, InterruptT>& parent)
: mParent(&parent)
, mMinAbsS(ValueType(1e-6))
, mMap(parent.mTracker.grid().transform().template constMap<MapT>().get())
, mTask(nullptr)
{
}
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
Morph(const Morph& other)
: mParent(other.mParent)
, mMinAbsS(other.mMinAbsS)
, mMaxAbsS(other.mMaxAbsS)
, mMap(other.mMap)
, mTask(other.mTask)
{
}
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
Morph(Morph& other, tbb::split)
: mParent(other.mParent)
, mMinAbsS(other.mMinAbsS)
, mMaxAbsS(other.mMaxAbsS)
, mMap(other.mMap)
, mTask(other.mTask)
{
}
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline size_t
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
advect(ValueType time0, ValueType time1)
{
namespace ph = std::placeholders;
// Make sure we have enough temporal auxiliary buffers for the time
// integration AS WELL AS an extra buffer with the speed function!
static const Index auxBuffers = 1 + (TemporalScheme == math::TVD_RK3 ? 2 : 1);
size_t countCFL = 0;
while (time0 < time1 && mParent->mTracker.checkInterrupter()) {
mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers);
const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers);
if ( math::isZero(dt) ) break;//V is essentially zero so terminate
OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time
switch(TemporalScheme) {
case math::TVD_RK1:
// Perform one explicit Euler step: t1 = t0 + dt
// Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2);
// Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
this->cook(PARALLEL_FOR, 1);
break;
case math::TVD_RK2:
// Perform one explicit Euler step: t1 = t0 + dt
// Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]|
mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2);
// Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
this->cook(PARALLEL_FOR, 1);
// Convex combine explict Euler step: t2 = t0 + dt
// Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|)
mTask = std::bind(&Morph::euler12, ph::_1, ph::_2, dt);
// Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1)
this->cook(PARALLEL_FOR, 1);
break;
case math::TVD_RK3:
// Perform one explicit Euler step: t1 = t0 + dt
// Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]|
mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/3);
// Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1)
this->cook(PARALLEL_FOR, 1);
// Convex combine explict Euler step: t2 = t0 + dt/2
// Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|)
mTask = std::bind(&Morph::euler34, ph::_1, ph::_2, dt);
// Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2)
this->cook(PARALLEL_FOR, 2);
// Convex combine explict Euler step: t3 = t0 + dt
// Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|)
mTask = std::bind(&Morph::euler13, ph::_1, ph::_2, dt);
// Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2)
this->cook(PARALLEL_FOR, 2);
break;
case math::UNKNOWN_TIS:
default:
OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!");
}//end of compile-time resolved switch
OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
time0 += dt;
++countCFL;
mParent->mTracker.leafs().removeAuxBuffers();
// Track the narrow band
mParent->mTracker.track();
}//end wile-loop over time
return countCFL;//number of CLF propagation steps
}
template<typename GridT, typename InterruptT>
template<typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline typename GridT::ValueType
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer)
{
namespace ph = std::placeholders;
mMaxAbsS = mMinAbsS;
const size_t leafCount = mParent->mTracker.leafs().leafCount();
if (leafCount==0 || time0 >= time1) return ValueType(0);
const math::Transform& xform = mParent->mTracker.grid().transform();
if (mParent->mTarget->transform() == xform &&
(mParent->mMask == nullptr || mParent->mMask->transform() == xform)) {
mTask = std::bind(&Morph::sampleAlignedSpeed, ph::_1, ph::_2, speedBuffer);
} else {
mTask = std::bind(&Morph::sampleXformedSpeed, ph::_1, ph::_2, speedBuffer);
}
this->cook(PARALLEL_REDUCE);
if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ValueType(0);//speed is essentially zero
static const ValueType CFL = (TemporalScheme == math::TVD_RK1 ? ValueType(0.3) :
TemporalScheme == math::TVD_RK2 ? ValueType(0.9) :
ValueType(1.0))/math::Sqrt(ValueType(3.0));
const ValueType dt = math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize();
return math::Min(dt, ValueType(CFL*dx/mMaxAbsS));
}
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline void
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
{
using VoxelIterT = typename LeafType::ValueOnCIter;
using SamplerT = tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler>;
const MapT& map = *mMap;
mParent->mTracker.checkInterrupter();
typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor();
SamplerT target(targetAcc, mParent->mTarget->transform());
if (mParent->mMask == nullptr) {
for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
ValueType* speed = leafIter.buffer(speedBuffer).data();
bool isZero = true;
for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
ValueType& s = speed[voxelIter.pos()];
s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
if (!math::isApproxZero(s)) isZero = false;
mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
}
if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
}
} else {
const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
const bool invMask = mParent->isMaskInverted();
typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor();
SamplerT mask(maskAcc, mParent->mMask->transform());
for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
ValueType* speed = leafIter.buffer(speedBuffer).data();
bool isZero = true;
for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());//world space
const ValueType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm);
ValueType& s = speed[voxelIter.pos()];
s -= target.wsSample(xyz);
s *= invMask ? 1 - a : a;
if (!math::isApproxZero(s)) isZero = false;
mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
}
if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
}
}
}
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline void
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
{
using VoxelIterT = typename LeafType::ValueOnCIter;
mParent->mTracker.checkInterrupter();
typename GridT::ConstAccessor target = mParent->mTarget->getAccessor();
if (mParent->mMask == nullptr) {
for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
ValueType* speed = leafIter.buffer(speedBuffer).data();
bool isZero = true;
for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
ValueType& s = speed[voxelIter.pos()];
s -= target.getValue(voxelIter.getCoord());
if (!math::isApproxZero(s)) isZero = false;
mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
}
if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
}
} else {
const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask);
const bool invMask = mParent->isMaskInverted();
typename GridT::ConstAccessor mask = mParent->mMask->getAccessor();
for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
ValueType* speed = leafIter.buffer(speedBuffer).data();
bool isZero = true;
for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
const Coord ijk = voxelIter.getCoord();//index space
const ValueType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm);
ValueType& s = speed[voxelIter.pos()];
s -= target.getValue(ijk);
s *= invMask ? 1 - a : a;
if (!math::isApproxZero(s)) isZero = false;
mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
}
if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
}
}
}
template<typename GridT, typename InterruptT>
template <typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
inline void
LevelSetMorphing<GridT, InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
cook(ThreadingMode mode, size_t swapBuffer)
{
mParent->mTracker.startInterrupter("Morphing level set");
const int grainSize = mParent->mTracker.getGrainSize();
const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize);
if (mParent->mTracker.getGrainSize()==0) {
(*this)(range);
} else if (mode == PARALLEL_FOR) {
tbb::parallel_for(range, *this);
} else if (mode == PARALLEL_REDUCE) {
tbb::parallel_reduce(range, *this);
} else {
OPENVDB_THROW(ValueError, "expected threading mode " << int(PARALLEL_FOR)
<< " or " << int(PARALLEL_REDUCE) << ", got " << int(mode));
}
mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0);
mParent->mTracker.endInterrupter();
}
template<typename GridT, typename InterruptT>
template<typename MapT, math::BiasedGradientScheme SpatialScheme,
math::TemporalIntegrationScheme TemporalScheme>
template <int Nominator, int Denominator>
inline void
LevelSetMorphing<GridT,InterruptT>::
Morph<MapT, SpatialScheme, TemporalScheme>::
euler(const LeafRange& range, ValueType dt,
Index phiBuffer, Index resultBuffer, Index speedBuffer)
{
using SchemeT = math::BIAS_SCHEME<SpatialScheme>;
using StencilT = typename SchemeT::template ISStencil<GridType>::StencilType;
using VoxelIterT = typename LeafType::ValueOnCIter;
using NumGrad = math::GradientNormSqrd<MapT, SpatialScheme>;
static const ValueType Alpha = ValueType(Nominator)/ValueType(Denominator);
static const ValueType Beta = ValueType(1) - Alpha;
mParent->mTracker.checkInterrupter();
const MapT& map = *mMap;
StencilT stencil(mParent->mTracker.grid());
for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
const ValueType* speed = leafIter.buffer(speedBuffer).data();
if (math::isExactlyEqual(speed[0], std::numeric_limits<ValueType>::max())) continue;
const ValueType* phi = leafIter.buffer(phiBuffer).data();
ValueType* result = leafIter.buffer(resultBuffer).data();
for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
const Index n = voxelIter.pos();
if (math::isApproxZero(speed[n])) continue;
stencil.moveTo(voxelIter);
const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
}//loop over active voxels in the leaf of the mask
}//loop over leafs of the level set
}
////////////////////////////////////////
// Explicit Template Instantiation
#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
#ifdef OPENVDB_INSTANTIATE_LEVELSETMORPH
#include <openvdb/util/ExplicitInstantiation.h>
#endif
OPENVDB_INSTANTIATE_CLASS LevelSetMorphing<FloatGrid, util::NullInterrupter>;
OPENVDB_INSTANTIATE_CLASS LevelSetMorphing<DoubleGrid, util::NullInterrupter>;
#endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
} // namespace tools
} // namespace OPENVDB_VERSION_NAME
} // namespace openvdb
#endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED
ExplicitInstantiation.h
FiniteDifference.h
Interpolation.h
LevelSetTracker.h
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition Platform.h:118
OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition Platform.h:117
openvdb::v10_0::ValueError
Definition Exceptions.h:65
openvdb::v10_0::math::Transform
Definition Transform.h:40
openvdb::v10_0::math::Transform::mapType
Name mapType() const
Return the transformation map's type-name.
Definition Transform.h:66
openvdb::v10_0::math::UniformScaleMap::mapType
static Name mapType()
Definition Maps.h:939
openvdb::v10_0::tools::LevelSetMorphing
Shape morphology of level sets. Morphing from a source narrow-band level sets to a target narrow-band...
Definition LevelSetMorph.h:48
openvdb::v10_0::tools::LevelSetMorphing::LeafRange
typename TrackerT::LeafRange LeafRange
Definition LevelSetMorph.h:53
openvdb::v10_0::tools::LevelSetMorphing::minMask
ValueType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value.
Definition LevelSetMorph.h:122
openvdb::v10_0::tools::LevelSetMorphing::setTrackerSpatialScheme
void setTrackerSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition LevelSetMorph.h:95
openvdb::v10_0::tools::LevelSetMorphing::getTrackerSpatialScheme
math::BiasedGradientScheme getTrackerSpatialScheme() const
Return the spatial finite-difference scheme.
Definition LevelSetMorph.h:90
openvdb::v10_0::tools::LevelSetMorphing::setTrackerTemporalScheme
void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition LevelSetMorph.h:105
openvdb::v10_0::tools::LevelSetMorphing::setMaskRange
void setMaskRange(ValueType min, ValueType max)
Define the range for the (optional) scalar mask.
Definition LevelSetMorph.h:135
openvdb::v10_0::tools::LevelSetMorphing::ValueType
typename TrackerT::ValueType ValueType
Definition LevelSetMorph.h:56
openvdb::v10_0::tools::LevelSetMorphing::getTemporalScheme
math::TemporalIntegrationScheme getTemporalScheme() const
Return the temporal integration scheme.
Definition LevelSetMorph.h:85
openvdb::v10_0::tools::LevelSetMorphing::BufferType
typename TrackerT::BufferType BufferType
Definition LevelSetMorph.h:55
openvdb::v10_0::tools::LevelSetMorphing::setTemporalScheme
void setTemporalScheme(math::TemporalIntegrationScheme scheme)
Set the temporal integration scheme.
Definition LevelSetMorph.h:87
openvdb::v10_0::tools::LevelSetMorphing::setSpatialScheme
void setSpatialScheme(math::BiasedGradientScheme scheme)
Set the spatial finite-difference scheme.
Definition LevelSetMorph.h:82
openvdb::v10_0::tools::LevelSetMorphing::getSpatialScheme
math::BiasedGradientScheme getSpatialScheme() const
Return the spatial finite-difference scheme.
Definition LevelSetMorph.h:80
openvdb::v10_0::tools::LevelSetMorphing::setGrainSize
void setGrainSize(int grainsize)
Set the grain size used for multithreading.
Definition LevelSetMorph.h:118
openvdb::v10_0::tools::LevelSetMorphing::invertMask
void invertMask(bool invert=true)
Invert the optional mask, i.e. min->max in the original mask maps to 1->0 in the inverted alpha mask.
Definition LevelSetMorph.h:147
openvdb::v10_0::tools::LevelSetMorphing::getGrainSize
int getGrainSize() const
Return the grain size used for multithreading.
Definition LevelSetMorph.h:115
openvdb::v10_0::tools::LevelSetMorphing::setNormCount
void setNormCount(int n)
Set the number of normalizations performed per track or normalize call.
Definition LevelSetMorph.h:112
openvdb::v10_0::tools::LevelSetMorphing::LeafType
typename TrackerT::LeafType LeafType
Definition LevelSetMorph.h:54
openvdb::v10_0::tools::LevelSetMorphing::maxMask
ValueType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value.
Definition LevelSetMorph.h:126
openvdb::v10_0::tools::LevelSetMorphing::TreeType
typename GridT::TreeType TreeType
Definition LevelSetMorph.h:51
openvdb::v10_0::tools::LevelSetMorphing::TrackerT
LevelSetTracker< GridT, InterruptT > TrackerT
Definition LevelSetMorph.h:52
openvdb::v10_0::tools::LevelSetMorphing::LevelSetMorphing
LevelSetMorphing(GridT &sourceGrid, const GridT &targetGrid, InterruptT *interrupt=nullptr)
Main constructor.
Definition LevelSetMorph.h:59
openvdb::v10_0::tools::LevelSetMorphing::getTrackerTemporalScheme
math::TemporalIntegrationScheme getTrackerTemporalScheme() const
Return the temporal integration scheme.
Definition LevelSetMorph.h:100
openvdb::v10_0::tools::LevelSetMorphing::setAlphaMask
void setAlphaMask(const GridT &maskGrid)
Define the alpha mask.
Definition LevelSetMorph.h:77
openvdb::v10_0::tools::LevelSetMorphing::advect
size_t advect(ValueType time0, ValueType time1)
Advect the level set from its current time, time0, to its final time, time1. If time0 > time1,...
Definition LevelSetMorph.h:241
openvdb::v10_0::tools::LevelSetMorphing::setTarget
void setTarget(const GridT &targetGrid)
Redefine the target level set.
Definition LevelSetMorph.h:74
openvdb::v10_0::tools::LevelSetMorphing::isMaskInverted
bool isMaskInverted() const
Return true if the mask is inverted, i.e. min->max in the original mask maps to 1->0 in the inverted ...
Definition LevelSetMorph.h:144
openvdb::v10_0::tools::LevelSetMorphing::GridType
GridT GridType
Definition LevelSetMorph.h:50
openvdb::v10_0::tools::LevelSetMorphing::getNormCount
int getNormCount() const
Return the number of normalizations performed per track or normalize call.
Definition LevelSetMorph.h:110
openvdb::v10_0::tools::LevelSetMorphing::~LevelSetMorphing
virtual ~LevelSetMorphing()
Definition LevelSetMorph.h:71
openvdb::v10_0::tools::LevelSetTracker
Performs multi-threaded interface tracking of narrow band level sets.
Definition LevelSetTracker.h:58
openvdb::v10_0::tools::LevelSetTracker::checkInterrupter
bool checkInterrupter()
Definition LevelSetTracker.h:385
openvdb::v10_0::tools::LevelSetTracker< GridT, InterruptT >::ValueType
typename TreeType::ValueType ValueType
Definition LevelSetTracker.h:65
openvdb::v10_0::tools::LevelSetTracker::startInterrupter
void startInterrupter(const char *msg)
Definition LevelSetTracker.h:369
openvdb::v10_0::tools::LevelSetTracker< GridT, InterruptT >::LeafRange
typename LeafManagerType::LeafRange LeafRange
Definition LevelSetTracker.h:67
openvdb::v10_0::tools::LevelSetTracker< GridT, InterruptT >::BufferType
typename LeafManagerType::BufferType BufferType
Definition LevelSetTracker.h:68
openvdb::v10_0::tools::LevelSetTracker::leafs
LeafManagerType & leafs()
Definition LevelSetTracker.h:189
openvdb::v10_0::tools::LevelSetTracker< GridT, InterruptT >::LeafType
typename TreeType::LeafNodeType LeafType
Definition LevelSetTracker.h:64
openvdb::v10_0::math
Definition Types.h:28
openvdb::v10_0::math::TemporalIntegrationScheme
TemporalIntegrationScheme
Temporal integration schemes.
Definition FiniteDifference.h:234
openvdb::v10_0::math::TVD_RK1
@ TVD_RK1
Definition FiniteDifference.h:236
openvdb::v10_0::math::UNKNOWN_TIS
@ UNKNOWN_TIS
Definition FiniteDifference.h:235
openvdb::v10_0::math::TVD_RK2
@ TVD_RK2
Definition FiniteDifference.h:237
openvdb::v10_0::math::TVD_RK3
@ TVD_RK3
Definition FiniteDifference.h:238
openvdb::v10_0::math::SmoothUnitStep
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else (3 − 2 x) x².
Definition Math.h:287
openvdb::v10_0::math::isApproxZero
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance.
Definition Math.h:349
openvdb::v10_0::math::Sqrt
float Sqrt(float x)
Return the square root of a floating-point value.
Definition Math.h:761
openvdb::v10_0::math::isApproxEqual
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition Math.h:406
openvdb::v10_0::math::Max
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition Math.h:595
openvdb::v10_0::math::BiasedGradientScheme
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition FiniteDifference.h:165
openvdb::v10_0::math::FIRST_BIAS
@ FIRST_BIAS
Definition FiniteDifference.h:167
openvdb::v10_0::math::THIRD_BIAS
@ THIRD_BIAS
Definition FiniteDifference.h:169
openvdb::v10_0::math::WENO5_BIAS
@ WENO5_BIAS
Definition FiniteDifference.h:170
openvdb::v10_0::math::UNKNOWN_BIAS
@ UNKNOWN_BIAS
Definition FiniteDifference.h:166
openvdb::v10_0::math::SECOND_BIAS
@ SECOND_BIAS
Definition FiniteDifference.h:168
openvdb::v10_0::math::HJWENO5_BIAS
@ HJWENO5_BIAS
Definition FiniteDifference.h:171
openvdb::v10_0::math::Min
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition Math.h:656
openvdb::v10_0::math::Abs
Coord Abs(const Coord &xyz)
Definition Coord.h:515
openvdb::v10_0::math::isExactlyEqual
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition Math.h:443
openvdb::v10_0::math::isZero
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition Math.h:337
openvdb::v10_0::tools::composite::min
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition Composite.h:106
openvdb::v10_0::Index
Index32 Index
Definition Types.h:54
openvdb::v10_0::Vec3R
math::Vec3< Real > Vec3R
Definition Types.h:72
openvdb
Definition Exceptions.h:13
OPENVDB_THROW
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
OPENVDB_VERSION_NAME
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
OPENVDB_USE_VERSION_NAMESPACE
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
OPENVDB_INSTANTIATE_CLASS
#define OPENVDB_INSTANTIATE_CLASS
Definition version.h.in:153