VolumeAdvect.h

Go to the documentation of this file.

// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0
//
//

#ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
#define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED

#include <tbb/parallel_for.h>
#include <openvdb/Types.h>
#include <openvdb/math/Math.h>
#include <openvdb/util/NullInterrupter.h>
#include "Interpolation.h"// for Sampler
#include "VelocityFields.h" // for VelocityIntegrator
#include "Morphology.h"//for dilateActiveValues and dilateVoxels
#include "Prune.h"// for prune
#include "Statistics.h" // for extrema
#include <functional>


namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace tools {

namespace Scheme {
    enum SemiLagrangian { SEMI, MID, RK3, RK4, MAC, BFECC };
    enum Limiter { NO_LIMITER, CLAMP, REVERT };
}


template<typename VelocityGridT = Vec3fGrid,
         bool StaggeredVelocity = false,
         typename InterrupterType = util::NullInterrupter>
class VolumeAdvection
{
public:

    VolumeAdvection(const VelocityGridT& velGrid, InterrupterType* interrupter = nullptr)
        : mVelGrid(velGrid)
        , mInterrupter(interrupter)
        , mIntegrator( Scheme::SEMI )
        , mLimiter( Scheme::CLAMP )
        , mGrainSize( 128 )
        , mSubSteps( 1 )
    {
        math::Extrema e = extrema(velGrid.cbeginValueAll(), /*threading*/true);
        e.add(velGrid.background().length());
        mMaxVelocity = e.max();
    }

    virtual ~VolumeAdvection()
    {
    }

    int spatialOrder() const { return (mIntegrator == Scheme::MAC ||
                                       mIntegrator == Scheme::BFECC) ? 2 : 1; }

    int temporalOrder() const {
        switch (mIntegrator) {
        case Scheme::SEMI: return 1;
        case Scheme::MID:  return 2;
        case Scheme::RK3:  return 3;
        case Scheme::RK4:  return 4;
        case Scheme::BFECC:return 2;
        case Scheme::MAC:  return 2;
        }
        return 0;//should never reach this point
    }

    void setIntegrator(Scheme::SemiLagrangian integrator) { mIntegrator = integrator; }

    Scheme::SemiLagrangian getIntegrator() const { return mIntegrator; }

    void setLimiter(Scheme::Limiter limiter) { mLimiter = limiter; }

    Scheme::Limiter getLimiter() const { return mLimiter; }

    bool isLimiterOn() const { return this->spatialOrder()>1 &&
                                      mLimiter != Scheme::NO_LIMITER; }

    size_t getGrainSize() const { return mGrainSize; }

    void setGrainSize(size_t grainsize) { mGrainSize = grainsize; }

    int getSubSteps() const { return mSubSteps; }

    void setSubSteps(int substeps) { mSubSteps = math::Max(1, substeps); }

    double getMaxVelocity() const { return mMaxVelocity; }

    template<typename VolumeGridT>
    int getMaxDistance(const VolumeGridT& inGrid, double dt) const
    {
        if (!inGrid.hasUniformVoxels()) {
            OPENVDB_THROW(RuntimeError, "Volume grid does not have uniform voxels!");
        }
        const double d = mMaxVelocity*math::Abs(dt)/inGrid.voxelSize()[0];
        return static_cast<int>( math::RoundUp(d) );
    }

    template<typename VolumeGridT,
             typename VolumeSamplerT>//only C++11 allows for a default argument
    typename VolumeGridT::Ptr advect(const VolumeGridT& inGrid, double timeStep)
    {
        typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
        const double dt = timeStep/mSubSteps;
        const int n = this->getMaxDistance(inGrid, dt);
        dilateActiveValues( outGrid->tree(), n, NN_FACE, EXPAND_TILES);
        this->template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
        for (int step = 1; step < mSubSteps; ++step) {
            typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
            dilateActiveValues( tmpGrid->tree(), n, NN_FACE, EXPAND_TILES);
            this->template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
            outGrid.swap( tmpGrid );
        }

        return outGrid;
    }

    template<typename VolumeGridT,
             typename MaskGridT,
             typename VolumeSamplerT>//only C++11 allows for a default argument
    typename VolumeGridT::Ptr advect(const VolumeGridT& inGrid, const MaskGridT& mask, double timeStep)
    {
        if (inGrid.transform() != mask.transform()) {
            OPENVDB_THROW(RuntimeError, "Volume grid and mask grid are misaligned! Consider "
                          "resampling either of the two grids into the index space of the other.");
        }
        typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
        const double dt = timeStep/mSubSteps;
        const int n = this->getMaxDistance(inGrid, dt);
        dilateActiveValues( outGrid->tree(), n, NN_FACE, EXPAND_TILES);
        outGrid->topologyIntersection( mask );
        pruneInactive( outGrid->tree(), mGrainSize>0, mGrainSize );
        this->template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
        outGrid->topologyUnion( inGrid );

        for (int step = 1; step < mSubSteps; ++step) {
            typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
            dilateActiveValues( tmpGrid->tree(), n, NN_FACE, EXPAND_TILES);
            tmpGrid->topologyIntersection( mask );
            pruneInactive( tmpGrid->tree(), mGrainSize>0, mGrainSize );
            this->template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
            tmpGrid->topologyUnion( inGrid );
            outGrid.swap( tmpGrid );
        }
        return outGrid;
    }

private:
    // disallow copy construction and copy by assignment!
    VolumeAdvection(const VolumeAdvection&);// not implemented
    VolumeAdvection& operator=(const VolumeAdvection&);// not implemented

    void start(const char* str) const
    {
        if (mInterrupter) mInterrupter->start(str);
    }
    void stop() const
    {
        if (mInterrupter) mInterrupter->end();
    }
    bool interrupt() const
    {
        if (mInterrupter && util::wasInterrupted(mInterrupter)) {
            tbb::task::self().cancel_group_execution();
            return true;
        }
        return false;
    }

    template<typename VolumeGridT, typename VolumeSamplerT>
    void cook(VolumeGridT& outGrid, const VolumeGridT& inGrid, double dt)
    {
        switch (mIntegrator) {
        case Scheme::SEMI: {
            Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *this);
            adv.cook(outGrid, dt);
            break;
        }
        case Scheme::MID: {
            Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *this);
            adv.cook(outGrid, dt);
            break;
        }
        case Scheme::RK3: {
            Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *this);
            adv.cook(outGrid, dt);
            break;
        }
        case Scheme::RK4: {
            Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *this);
            adv.cook(outGrid, dt);
            break;
        }
        case Scheme::BFECC: {
            Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *this);
            adv.cook(outGrid, dt);
            break;
        }
        case Scheme::MAC: {
            Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *this);
            adv.cook(outGrid, dt);
            break;
        }
        default:
            OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!");
        }
        pruneInactive(outGrid.tree(), mGrainSize>0, mGrainSize);
    }

    // Private class that implements the multi-threaded advection
    template<typename VolumeGridT, size_t OrderRK, typename SamplerT> struct Advect;

    // Private member data of VolumeAdvection
    const VelocityGridT&   mVelGrid;
    double                 mMaxVelocity;
    InterrupterType*       mInterrupter;
    Scheme::SemiLagrangian mIntegrator;
    Scheme::Limiter        mLimiter;
    size_t                 mGrainSize;
    int                    mSubSteps;
};//end of VolumeAdvection class

// Private class that implements the multi-threaded advection
template<typename VelocityGridT, bool StaggeredVelocity, typename InterrupterType>
template<typename VolumeGridT, size_t OrderRK, typename SamplerT>
struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
{
    using TreeT = typename VolumeGridT::TreeType;
    using AccT = typename VolumeGridT::ConstAccessor;
    using ValueT = typename TreeT::ValueType;
    using LeafManagerT = typename tree::LeafManager<TreeT>;
    using LeafNodeT = typename LeafManagerT::LeafNodeType;
    using LeafRangeT = typename LeafManagerT::LeafRange;
    using VelocityIntegratorT = VelocityIntegrator<VelocityGridT, StaggeredVelocity>;
    using RealT = typename VelocityIntegratorT::ElementType;
    using VoxelIterT = typename TreeT::LeafNodeType::ValueOnIter;

    Advect(const VolumeGridT& inGrid, const VolumeAdvection& parent)
        : mTask(nullptr)
        , mInGrid(&inGrid)
        , mVelocityInt(parent.mVelGrid)
        , mParent(&parent)
    {
    }
    inline void cook(const LeafRangeT& range)
    {
        if (mParent->mGrainSize > 0) {
            tbb::parallel_for(range, *this);
        } else {
            (*this)(range);
        }
    }
    void operator()(const LeafRangeT& range) const
    {
        assert(mTask);
        mTask(const_cast<Advect*>(this), range);
    }
    void cook(VolumeGridT& outGrid, double time_step)
    {
        namespace ph = std::placeholders;

        mParent->start("Advecting volume");
        LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
        const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
        const RealT dt = static_cast<RealT>(-time_step);//method of characteristics backtracks
        if (mParent->mIntegrator == Scheme::MAC) {
            mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);//out[0]=forward
            this->cook(range);
            mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);//out[1]=backward
            this->cook(range);
            mTask = std::bind(&Advect::mac, ph::_1, ph::_2);//out[0] = out[0] + (in[0] - out[1])/2
            this->cook(range);
        } else if (mParent->mIntegrator == Scheme::BFECC) {
            mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);//out[0]=forward
            this->cook(range);
            mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);//out[1]=backward
            this->cook(range);
            mTask = std::bind(&Advect::bfecc, ph::_1, ph::_2);//out[0] = (3*in[0] - out[1])/2
            this->cook(range);
            mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 1, &outGrid);//out[1]=forward
            this->cook(range);
            manager.swapLeafBuffer(1);// out[0] = out[1]
        } else {// SEMI, MID, RK3 and RK4
            mTask = std::bind(&Advect::rk, ph::_1, ph::_2,  dt, 0, mInGrid);//forward
            this->cook(range);
        }

        if (mParent->spatialOrder()==2) manager.removeAuxBuffers();

        mTask = std::bind(&Advect::limiter, ph::_1, ph::_2, dt);// out[0] = limiter( out[0] )
        this->cook(range);

        mParent->stop();
    }
    // Last step of the MacCormack scheme: out[0] = out[0] + (in[0] - out[1])/2
    void mac(const LeafRangeT& range) const
    {
        if (mParent->interrupt()) return;
        assert( mParent->mIntegrator == Scheme::MAC );
        AccT acc = mInGrid->getAccessor();
        for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
            ValueT* out0 = leafIter.buffer( 0 ).data();// forward
            const ValueT* out1 = leafIter.buffer( 1 ).data();// backward
            const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
            if (leaf != nullptr) {
                const ValueT* in0 = leaf->buffer().data();
                for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
                    const Index i = voxelIter.pos();
                    out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
                }
            } else {
                for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
                    const Index i = voxelIter.pos();
                    out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
                }//loop over active voxels
            }
        }//loop over leaf nodes
    }
    // Intermediate step in the BFECC scheme: out[0] = (3*in[0] - out[1])/2
    void bfecc(const LeafRangeT& range) const
    {
        if (mParent->interrupt()) return;
        assert( mParent->mIntegrator == Scheme::BFECC );
        AccT acc = mInGrid->getAccessor();
        for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
            ValueT* out0 = leafIter.buffer( 0 ).data();// forward
            const ValueT* out1 = leafIter.buffer( 1 ).data();// backward
            const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
            if (leaf != nullptr) {
                const ValueT* in0 = leaf->buffer().data();
                for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
                    const Index i = voxelIter.pos();
                    out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
                }//loop over active voxels
            } else {
                for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
                    const Index i = voxelIter.pos();
                    out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
                }//loop over active voxels
            }
        }//loop over leaf nodes
    }
    // Semi-Lagrangian integration with Runge-Kutta of various orders (1->4)
    void rk(const LeafRangeT& range, RealT dt, size_t n, const VolumeGridT* grid) const
    {
        if (mParent->interrupt()) return;
        const math::Transform& xform = mInGrid->transform();
        AccT acc = grid->getAccessor();
        for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
            ValueT* phi = leafIter.buffer( n ).data();
            for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
                ValueT& value = phi[voxelIter.pos()];
                Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
                mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
                value = SamplerT::sample(acc, xform.worldToIndex(wPos));
            }//loop over active voxels
        }//loop over leaf nodes
    }
    void limiter(const LeafRangeT& range, RealT dt) const
    {
        if (mParent->interrupt()) return;
        const bool doLimiter = mParent->isLimiterOn();
        const bool doClamp = mParent->mLimiter == Scheme::CLAMP;
        ValueT data[2][2][2], vMin, vMax;
        const math::Transform& xform = mInGrid->transform();
        AccT acc = mInGrid->getAccessor();
        const ValueT backg = mInGrid->background();
        for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
            ValueT* phi = leafIter.buffer( 0 ).data();
            for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
                ValueT& value = phi[voxelIter.pos()];

                if ( doLimiter ) {
                    assert(OrderRK == 1);
                    Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
                    mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);// Explicit Euler
                    Vec3d iPos = xform.worldToIndex(wPos);
                    Coord ijk  = Coord::floor( iPos );
                    BoxSampler::getValues(data, acc, ijk);
                    BoxSampler::extrema(data, vMin, vMax);
                    if ( doClamp ) {
                        value = math::Clamp( value, vMin, vMax);
                    } else if (value < vMin || value > vMax ) {
                        iPos -= Vec3R(ijk[0], ijk[1], ijk[2]);//unit coordinates
                        value = BoxSampler::trilinearInterpolation( data, iPos );
                    }
                }

                if (math::isApproxEqual(value, backg, math::Delta<ValueT>::value())) {
                    value = backg;
                    leafIter->setValueOff( voxelIter.pos() );
                }
            }//loop over active voxels
        }//loop over leaf nodes
    }
    // Public member data of the private Advect class

    typename std::function<void (Advect*, const LeafRangeT&)> mTask;
    const VolumeGridT*        mInGrid;
    const VelocityIntegratorT mVelocityInt;// lightweight!
    const VolumeAdvection*    mParent;
};// end of private member class Advect

} // namespace tools
} // namespace OPENVDB_VERSION_NAME
} // namespace openvdb

#endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED