fastmarch.h

Go to the documentation of this file.

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey 
 *
 * This program is free software, distributed under the terms of the
 * GNU General Public License (GPL) 
 * http://www.gnu.org/licenses
 *
 * Fast marching
 *
 ******************************************************************************/

#ifndef _FASTMARCH_H
#define _FASTMARCH_H

#include <queue>
#include "levelset.h"

namespace Manta {
    
//  This class exists in two versions: for scalar, and for vector values - the only difference are 
//  flag checks i transpTouch (for simplicity in separate classes)

template<class GRID, class T>
inline T fmInterpolateNeighbors(GRID* mpVal, int x,int y,int z, Real *weights) {
    T val(0.); 
    if(weights[0]>0.0) val += mpVal->get(x+1, y+0, z+0) * weights[0];
    if(weights[1]>0.0) val += mpVal->get(x-1, y+0, z+0) * weights[1];
    if(weights[2]>0.0) val += mpVal->get(x+0, y+1, z+0) * weights[2];
    if(weights[3]>0.0) val += mpVal->get(x+0, y-1, z+0) * weights[3];
    if(mpVal->is3D()) {
        if(weights[4]>0.0) val += mpVal->get(x+0, y+0, z+1) * weights[4];
        if(weights[5]>0.0) val += mpVal->get(x+0, y+0, z-1) * weights[5];
    } 
    return val;
}

template<class GRID, class T>
class FmValueTransportScalar {
public:
    FmValueTransportScalar() : mpVal(0),mpFlags(0) { };
    ~FmValueTransportScalar() { }; 
    void initMarching(GRID* val, FlagGrid* flags) {
        mpVal = val;
        mpFlags = flags;
    } 
    inline bool isInitialized() { return mpVal != 0; }

    inline void transpTouch(int x,int y,int z, Real *weights, Real time) {
        if(!mpVal || !mpFlags->isEmpty(x,y,z)) return;
        T val = fmInterpolateNeighbors<GRID,T>(mpVal,x,y,z,weights); 
        (*mpVal)(x,y,z) = val;
    }; 
protected:
    GRID* mpVal;
    FlagGrid* mpFlags;
};

template<class GRID, class T>
class FmValueTransportVec3 {
public:
    FmValueTransportVec3() : mpVal(0), mpFlags(0) { };
    ~FmValueTransportVec3() { };
    inline bool isInitialized() { return mpVal != 0; } 
    void initMarching(GRID* val, const FlagGrid* flags) {
        mpVal = val;
        mpFlags = flags;
    } 

    inline void transpTouch(int x,int y,int z, Real *weights, Real time) {
        if(!mpVal || !mpFlags->isEmpty(x,y,z)) return;
        
        T val = fmInterpolateNeighbors<GRID,T>(mpVal,x,y,z,weights); 

        // set velocity components if adjacent is empty
        if (mpFlags->isEmpty(x-1,y,z)) (*mpVal)(x,y,z).x = val.x;
        if (mpFlags->isEmpty(x,y-1,z)) (*mpVal)(x,y,z).y = val.y;
        if(mpVal->is3D()) { if (mpFlags->isEmpty(x,y,z-1)) (*mpVal)(x,y,z).z = val.z; }
    }; 

protected:
    GRID* mpVal;
    const FlagGrid* mpFlags;
};

class FmHeapEntryOut {
public:
    Vec3i p;
    // quick time access for sorting
    Real time;
    static inline bool compare(const Real x, const Real y) { 
        return x > y;
    }

    inline bool operator< (const FmHeapEntryOut& o) const {
        const Real d = fabs((time) - ((o.time)));
        if (d > 0.) return (time) > ((o.time)); 
        if (p.z != o.p.z) return p.z > o.p.z;
        if (p.y != o.p.y) return p.y > o.p.y;
        return p.x > o.p.x;
    };

};

class FmHeapEntryIn {
public:
    Vec3i p;
    // quick time access for sorting
    Real time;
    static inline bool compare(const Real x, const Real y) { 
        return x < y;
    }

    inline bool operator< (const FmHeapEntryIn& o) const {
        const Real d = fabs((time) - ((o.time)));
        if (d > 0.) return (time) < ((o.time)); 
        if (p.z != o.p.z) return p.z < o.p.z;
        if (p.y != o.p.y) return p.y < o.p.y;
        return p.x < o.p.x;
    };
};


template<class T, int TDIR>
class FastMarch {

public:
    // MSVC doesn't allow static const variables in template classes
    static inline Real InvalidTime() { return -1000; }
    static inline Real InvtOffset() { return 500; }

    enum SpecialValues { FlagInited = 1, FlagIsOnHeap = 2};

    FastMarch(const FlagGrid& flags, Grid<int>& fmFlags, Grid<Real>& levelset, Real maxTime, MACGrid* velTransport = NULL);
    ~FastMarch() {}
    
    void performMarching();

    inline bool isInvalid(Real v) const { return (v <= InvalidTime()); }

    void addToList(const Vec3i& p, const Vec3i& src);

    inline Real phi2time(Real phival) { return (phival-InvalidTime()+ InvtOffset()) * -1.0; }
    
    inline Real time2phi(Real tval) { return (InvalidTime() - InvtOffset() - tval); }

    inline Real _phi(int i, int j, int k) { return mLevelset(i,j,k); }
protected:   
    Grid<Real>&   mLevelset;
    const FlagGrid&     mFlags;
    Grid<int>&    mFmFlags;
    
    FmValueTransportVec3<MACGrid     , Vec3> mVelTransport;
    
    Real mMaxTime;

    std::priority_queue<T, std::vector<T>, std::less<T> > mHeap;
    Real mReheapVal;

    Real mWeights[6];

    template<int C> inline Real calcWeights(int& okCnt, int& invcnt, Real* v, const Vec3i& idx);
    
    inline Real calculateDistance(const Vec3i& pos);
};

} // namespace
#endif