00001
00009 #ifndef FLUXES_H
00010 #define FLUXES_H
00011
00012 #include <Globals.h>
00013 #include <Schemes.h>
00014 #include <Approximators.h>
00015
00016 namespace levelset {
00017
00018
00019
00020
00025 template<unsigned DIM, class Approximator, typename Pixel = float>
00026 class Fluxes
00027 {
00028 public:
00029 typedef BaseImage<DIM, Pixel> ImageType;
00030 typedef typename ImageType::Index Index;
00031 typedef Schemes<DIM,Pixel> Schemes;
00032
00035
00037 static Pixel reinitialization(const ImageType & phi, const Index & ind);
00039 static Pixel extension(const ImageType & phi, const ImageType & u, const Index & ind);
00041 static Pixel advection(const ImageType & u, const Index & ind, const Pixel v[]);
00042 #if 0
00044 Pixel grass_fire(const ImageType & phi, const Pixel beta, const Index & ind);
00046 Pixel anisotropic_eikonal(const ImageType & phi, const Pixel beta, const ImageType & D, const Index & ind);
00047 #endif
00048
00050 };
00051
00052
00054
00055
00056
00057 template <unsigned DIM, class Approximator, typename Pixel>
00058 Pixel Fluxes<DIM, Approximator, Pixel>::reinitialization(const ImageType & phi, const Index & ind)
00059 {
00060 const Pixel p = phi(ind);
00061 const Pixel s = p / sqrt( p*p + Globals::DELTA*Globals::DELTA );
00062 Pixel delta = 0;
00063
00064 if (s > 0) {
00065 for (int i=1 ; i<=DIM ; i++){
00066 const Pixel pf = Approximator::forward(phi,ind,i);
00067 const Pixel pb = Approximator::backward(phi,ind,i);
00068 if (pf < 0) delta += pf*pf;
00069 if (pb > 0) delta += pb*pb;
00070 }
00071 }
00072 else {
00073 for (int i=1 ; i<=DIM ; i++){
00074 const Pixel pf = Approximator::forward(phi,ind,i);
00075 const Pixel pb = Approximator::backward(phi,ind,i);
00076 if (pf > 0) delta += pf*pf;
00077 if (pb < 0) delta += pb*pb;
00078 }
00079 }
00080
00081 return - s * (sqrt(delta) - 1);
00082 }
00083
00084 template <unsigned DIM, class Approximator, typename Pixel>
00085 Pixel Fluxes<DIM, Approximator, Pixel>::extension(const ImageType & phi, const ImageType & u, const Index & ind)
00086 {
00087 Pixel norm[DIM];
00088 Schemes::normal(phi, ind, norm);
00089 const Pixel s = (phi(ind) <= 0)?-1:1;
00090
00091 for (int i=0 ; i<DIM ; i++)
00092 norm[i] *= s;
00093
00094 return Fluxes<DIM, Approximator, Pixel>::advection(u, ind, norm);
00095 }
00096
00097 template <unsigned DIM, class Approximator, typename Pixel>
00098 Pixel Fluxes<DIM, Approximator, Pixel>::advection(const ImageType & u, const Index & ind, const Pixel v[])
00099 {
00100 Pixel result = 0;
00101
00102 for (int i=1 ; i<=DIM ; i++)
00103 result += v[i-1] * ( (v[i-1]>0) ? Approximator::backward(u,ind,i) : Approximator::forward(u,ind,i) );
00104
00105 return -result;
00106 }
00107
00108
00109 #if 0
00111 template <class approximator, typename Pixel>
00112 Pixel GrassFire(const ImageType & phi, const Pixel beta, const int x, const int y, const int z=0)
00113 {
00114
00115 const Pixel upx = approximator::upx(u,x,y);
00116 const Pixel umx = approximator::umx(u,x,y);
00117 const Pixel upy = approximator::upy(u,x,y);
00118 const Pixel umy = approximator::umy(u,x,y);
00119
00120 Pixel delta = 0;
00121
00122 if (beta > 0) {
00123 if (upx < 0) delta += upx*upx;
00124 if (umx > 0) delta += umx*umx;
00125 if (upy < 0) delta += upy*upy;
00126 if (umy > 0) delta += umy*umy;
00127 }
00128 else {
00129 if (upx > 0) delta += upx*upx;
00130 if (umx < 0) delta += umx*umx;
00131 if (upy > 0) delta += upy*upy;
00132 if (umy < 0) delta += umy*umy;
00133 }
00134
00135 return - beta * sqrt(delta);
00136
00137
00138 const Pixel upx = approximator::upx(u,x,y,z);
00139 const Pixel umx = approximator::umx(u,x,y,z);
00140 const Pixel upy = approximator::upy(u,x,y,z);
00141 const Pixel umy = approximator::umy(u,x,y,z);
00142 const Pixel upz = approximator::upz(u,x,y,z);
00143 const Pixel umz = approximator::umz(u,x,y,z);
00144
00145 Pixel delta = 0;
00146
00147 if (beta > 0) {
00148 if (upx < 0) delta += upx*upx;
00149 if (umx > 0) delta += umx*umx;
00150 if (upy < 0) delta += upy*upy;
00151 if (umy > 0) delta += umy*umy;
00152 if (upz < 0) delta += upz*upz;
00153 if (umz > 0) delta += umz*umz;
00154 }
00155 else {
00156 if (upx > 0) delta += upx*upx;
00157 if (umx < 0) delta += umx*umx;
00158 if (upy > 0) delta += upy*upy;
00159 if (umy < 0) delta += umy*umy;
00160 if (upz > 0) delta += upz*upz;
00161 if (umz < 0) delta += umz*umz;
00162 }
00163
00164 return - beta * sqrt(delta);
00165 }
00166
00167
00169 template <class approximator, typename Pixel>
00170 Pixel AnisotropicEikonal(const ImageType & phi, const Pixel beta, const MImage<Pixel, 2>& D, const int x, const int y, const int z=0)
00171 {
00172 const Pixel upx = approximator::upx(u,x,y);
00173 const Pixel umx = approximator::umx(u,x,y);
00174 const Pixel upy = approximator::upy(u,x,y);
00175 const Pixel umy = approximator::umy(u,x,y);
00176
00177 Pixel delta = 0.0;
00178 if (upx < 0) delta += D(0,0)*upx*upx;
00179 if (umx > 0) delta += D(0,0)*umx*umx;
00180 if (upy < 0) delta += D(1,1)*upy*upy;
00181 if (umy > 0) delta += D(1,1)*umy*umy;
00182
00183 delta += 2.0*D(0,1)*cimg::minmod(upx,umx)*cimg::minmod(upy,umy);
00184 return -beta*sqrt(delta);
00185
00186 const Pixel upx = approximator::upx(u,x,y,z);
00187 const Pixel umx = approximator::umx(u,x,y,z);
00188 const Pixel upy = approximator::upy(u,x,y,z);
00189 const Pixel umy = approximator::umy(u,x,y,z);
00190 const Pixel upz = approximator::upz(u,x,y,z);
00191 const Pixel umz = approximator::umz(u,x,y,z);
00192
00193 Pixel delta = 0.0;
00194 if (upx < 0) delta += D(0,0)*upx*upx;
00195 if (umx > 0) delta += D(0,0)*umx*umx;
00196 if (upy < 0) delta += D(1,1)*upy*upy;
00197 if (umy > 0) delta += D(1,1)*umy*umy;
00198 if (upz < 0) delta += D(2,2)*upz*upz;
00199 if (umz >0) delta += D(2,2)*umz*umz;
00200
00201 delta += 2.0*D(0,1)*minmod(upx,umx)*minmod(upy,umy);
00202 delta += 2.0*D(0,2)*minmod(upx,umx)*minmod(upz,umz);
00203 delta += 2.0*D(1,2)*minmod(upy,umy)*minmod(upz,umz);
00204 return -beta*sqrt(delta);
00205 }
00206
00207 #endif
00208
00209 }
00210
00211 #endif // FLUXES_H