text: /******************************************************************************
   *
   * MantaFlow fluid solver framework
   * Copyright 2011 Tobias Pfaff, Nils Thuerey 
   *
   * This program is free software, distributed under the terms of the
   * Apache License, Version 2.0 
   * http://www.apache.org/licenses/LICENSE-2.0
   *
   * Fast marching
   *
   ******************************************************************************/
  
  #ifndef _FASTMARCH_H
  #define _FASTMARCH_H
  
  #include <queue>
  #include "levelset.h"
  
  namespace Manta {
  	
  //! Fast marching. Transport certain values
  //  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; }
  
  	//! cell is touched by marching from source cell
  	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;
  	} 
  
  	//! cell is touched by marching from source cell
  	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;
  	};
  };
  
  
  //! fast marching algorithm wrapper class
  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() {}
  	
  	//! advect level set function with given velocity */
  	void performMarching();
  
  	//! test value for invalidity
  	inline bool isInvalid(Real v) const { return (v <= InvalidTime()); }
  
  	void addToList(const Vec3i& p, const Vec3i& src);
  
  	//! convert phi to time value
  	inline Real phi2time(Real phival) { return (phival-InvalidTime()+ InvtOffset()) * -1.0; }
  	
  	//! ... and back
  	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;
  	
  	//! velocity extrpolation
  	FmValueTransportVec3<MACGrid     , Vec3> mVelTransport;
  	
  	//! maximal time to march for
  	Real mMaxTime;
  
  	//! fast marching list
  	std::priority_queue<T, std::vector<T>, std::less<T> > mHeap;
  	Real mReheapVal;
  
  	//! weights for touching points
  	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
  
  