00001
00009 #ifndef APPROXIMATORS_H
00010 #define APPROXIMATORS_H
00011
00012 namespace levelset {
00013
00014
00015
00020 template<unsigned DIM, typename Pixel = float>
00021 class Weno3
00022 {
00023 public:
00024 typedef BaseImage<DIM, Pixel> ImageType;
00025 typedef typename ImageType::Index Index;
00026
00028 static Pixel forward(const ImageType & phi, const Index & ind, unsigned axis);
00030 static Pixel backward(const ImageType & phi, const Index & ind, unsigned axis);
00031 private:
00033 static Pixel combine(const Pixel v1, const Pixel v2, const Pixel v3);
00034 };
00035
00036
00038
00039
00040
00041 template<unsigned DIM, typename Pixel>
00042 Pixel Weno3<DIM, Pixel>::forward(const ImageType & phi, const Index & ind, unsigned axis)
00043 {
00044 int pos = ind(axis);
00045 int limit = phi.extent(axis-1)-1;
00046 Index ind_tmp = ind;
00047
00048 ind_tmp(axis) = (pos+2 > limit) ? limit : pos+2;
00049 const Pixel u3 = phi(ind_tmp);
00050
00051 ind_tmp(axis) = (pos+1 > limit) ? limit : pos+1;
00052 const Pixel u2 = phi(ind_tmp);
00053
00054 ind_tmp(axis) = pos;
00055 const Pixel u1 = phi(ind_tmp);
00056
00057 ind_tmp(axis) = (pos-1 < 0) ? 0 : pos-1;
00058 const Pixel u0 = phi(ind_tmp);
00059
00060 return combine(u3-u2,u2-u1,u1-u0);
00061 }
00062
00063 template<unsigned DIM, typename Pixel>
00064 Pixel Weno3<DIM, Pixel>::backward(const ImageType & phi, const Index & ind, unsigned axis)
00065 {
00066 int pos = ind(axis);
00067 int limit = phi.extent(axis-1)-1;
00068 Index ind_tmp = ind;
00069
00070 ind_tmp(axis) = (pos+1 > limit) ? limit : pos+1;
00071 const Pixel u3 = phi(ind_tmp);
00072
00073 ind_tmp(axis) = pos;
00074 const Pixel u2 = phi(ind_tmp);
00075
00076 ind_tmp(axis) = (pos-1 < 0) ? 0 : pos-1;
00077 const Pixel u1 = phi(ind_tmp);
00078
00079 ind_tmp(axis) = (pos-2 < 0) ? 0 : pos-2;
00080 const Pixel u0 = phi(ind_tmp);
00081
00082 return combine(u1-u0,u2-u1,u3-u2);
00083 }
00084
00085 template<unsigned DIM, typename Pixel>
00086 Pixel Weno3<DIM, Pixel>::combine(const Pixel v1, const Pixel v2, const Pixel v3)
00087 {
00088 Pixel s = ( Globals::EPS + Maths::sqr(v2 - v1) ) / ( Globals::EPS + Maths::sqr(v3 - v2) );
00089 s = 1 / (1 + 2*s*s);
00090 return (v2+v3 - s * (v1 - 2*v2 + v3)) / 2;
00091 }
00092
00093 #if 0
00094
00099 class Centered2
00100 {
00101 public:
00102
00103 template<typename Pixel>
00104 static Pixel ux(const ImageType & phi, int x, int y = 0)
00105 {
00106 return (u.neumann_pix2d(x+1,y) - u.neumann_pix2d(x-1,y)) / 2;
00107 }
00108
00109 template<typename Pixel>
00110 static Pixel ux(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00111 {
00112 return (u.neumann_pix3d(x+1,y,z) - u.neumann_pix3d(x-1,y,z)) / 2;
00113 }
00114
00115 template<typename Pixel>
00116 static Pixel uy(const ImageType & phi, int x, int y)
00117 {
00118 return (u.neumann_pix2d(x,y+1) - u.neumann_pix2d(x,y-1)) / 2;
00119 }
00120
00121 template<typename Pixel>
00122 static Pixel uy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00123 {
00124 return (u.neumann_pix3d(x,y+1,z) - u.neumann_pix3d(x,y-1,z)) / 2;
00125 }
00126
00127 template<typename Pixel>
00128 static Pixel uz(const MImage<Pixel, 3> &u, int x, int y, int z)
00129 {
00130 return (u.neumann_pix3d(x,y,z+1) - u.neumann_pix3d(x,y,z-1)) / 2;
00131 }
00132
00133 template<typename Pixel>
00134 static Pixel upx(const ImageType & phi, int x, int y = 0)
00135 {
00136 return ux(u,x,y);
00137 }
00138
00139 template<typename Pixel>
00140 static Pixel upx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00141 {
00142 return ux(u,x,y,z);
00143 }
00144
00145 template<typename Pixel>
00146 static Pixel umx(const ImageType & phi, int x, int y = 0)
00147 {
00148 return ux(u,x,y);
00149 }
00150
00151 template<typename Pixel>
00152 static Pixel umx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00153 {
00154 return ux(u,x,y,z);
00155 }
00156
00157 template<typename Pixel>
00158 static Pixel upy(const ImageType & phi, int x, int y)
00159 {
00160 return uy(u,x,y);
00161 }
00162
00163 template<typename Pixel>
00164 static Pixel upy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00165 {
00166 return uy(u,x,y,z);
00167 }
00168
00169 template<typename Pixel>
00170 static Pixel umy(const ImageType & phi, int x, int y)
00171 {
00172 return uy(u,x,y);
00173 }
00174
00175 template<typename Pixel>
00176 static Pixel umy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00177 {
00178 return uy(u,x,y,z);
00179 }
00180
00181 template<typename Pixel>
00182 static Pixel upz(const MImage<Pixel, 3> &u, int x, int y, int z)
00183 {
00184 return uz(u,x,y,z);
00185 }
00186
00187 template<typename Pixel>
00188 static Pixel umz(const MImage<Pixel, 3> &u, int x, int y, int z)
00189 {
00190 return uz(u,x,y,z);
00191 }
00192 };
00193
00194
00199 class Upwind
00200 {
00201 public:
00202
00203 template<typename Pixel>
00204 static Pixel upx(const ImageType & phi, int x, int y = 0)
00205 {
00206 return u.neumann_pix2d(x+1,y) - u.neumann_pix2d(x,y);
00207 }
00208
00209 template<typename Pixel>
00210 static Pixel upx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00211 {
00212 return u.neumann_pix3d(x+1,y,z) - u.neumann_pix3d(x,y,z);
00213 }
00214
00215 template<typename Pixel>
00216 static Pixel umx(const ImageType & phi, int x, int y = 0)
00217 {
00218 return u.neumann_pix2d(x,y) - u.neumann_pix2d(x-1,y);
00219 }
00220
00221 template<typename Pixel>
00222 static Pixel umx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00223 {
00224 return u.neumann_pix3d(x,y,z) - u.neumann_pix3d(x-1,y,z);
00225 }
00226
00227 template<typename Pixel>
00228 static Pixel upy(const ImageType & phi, int x, int y)
00229 {
00230 return u.neumann_pix2d(x,y+1) - u.neumann_pix2d(x,y);
00231 }
00232
00233 template<typename Pixel>
00234 static Pixel upy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00235 {
00236 return u.neumann_pix3d(x,y+1,z) - u.neumann_pix3d(x,y,z);
00237 }
00238
00239 template<typename Pixel>
00240 static Pixel umy(const ImageType & phi, int x, int y)
00241 {
00242 return u.neumann_pix2d(x,y) - u.neumann_pix2d(x,y-1);
00243 }
00244
00245 template<typename Pixel>
00246 static Pixel umy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00247 {
00248 return u.neumann_pix3d(x,y,z) - u.neumann_pix3d(x,y-1,z);
00249 }
00250
00251 template<typename Pixel>
00252 static Pixel upz(const MImage<Pixel, 3> &u, int x, int y, int z)
00253 {
00254 return u.neumann_pix3d(x,y,z+1) - u.neumann_pix3d(x,y,z);
00255 }
00256
00257 template<typename Pixel>
00258 static Pixel umz(const MImage<Pixel, 3> &u, int x, int y, int z)
00259 {
00260 return u.neumann_pix3d(x,y,z) - u.neumann_pix3d(x,y,z-1);
00261 }
00262 };
00263
00264
00269 class Weno3
00270 {
00271 public:
00272
00273 template<typename Pixel>
00274 static Pixel combine(const Pixel v1, const Pixel v2, const Pixel v3)
00275 {
00276 Pixel s = ( Globals::EPS + sqr(v2 - v1) ) / ( Globals::EPS + sqr(v3 - v2) );
00277 s = 1 / (1 + 2*s*s);
00278 return (v2+v3 - s * (v1 - 2*v2 + v3)) / 2;
00279 }
00280
00281 template<typename Pixel>
00282 static Pixel upx(const ImageType & phi, int x, int y = 0)
00283 {
00284 const Pixel u3 = u.neumann_pix2d(x+2,y);
00285 const Pixel u2 = u.neumann_pix2d(x+1,y);
00286 const Pixel u1 = u.neumann_pix2d(x,y);
00287 const Pixel u0 = u.neumann_pix2d(x-1,y);
00288 return combine(u3-u2,u2-u1,u1-u0);
00289 }
00290
00291 template<typename Pixel>
00292 static Pixel upx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00293 {
00294 const Pixel u3 = u.neumann_pix3d(x+2,y,z);
00295 const Pixel u2 = u.neumann_pix3d(x+1,y,z);
00296 const Pixel u1 = u.neumann_pix3d(x,y,z);
00297 const Pixel u0 = u.neumann_pix3d(x-1,y,z);
00298 return combine(u3-u2,u2-u1,u1-u0);
00299 }
00300
00301 template<typename Pixel>
00302 static Pixel umx(const ImageType & phi, int x, int y = 0)
00303 {
00304 const Pixel u3 = u.neumann_pix2d(x+1,y);
00305 const Pixel u2 = u.neumann_pix2d(x,y);
00306 const Pixel u1 = u.neumann_pix2d(x-1,y);
00307 const Pixel u0 = u.neumann_pix2d(x-2,y);
00308 return combine(u1-u0,u2-u1,u3-u2);
00309 }
00310
00311 template<typename Pixel>
00312 static Pixel umx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00313 {
00314 const Pixel u3 = u.neumann_pix3d(x+1,y,z);
00315 const Pixel u2 = u.neumann_pix3d(x,y,z);
00316 const Pixel u1 = u.neumann_pix3d(x-1,y,z);
00317 const Pixel u0 = u.neumann_pix3d(x-2,y,z);
00318 return combine(u1-u0,u2-u1,u3-u2);
00319 }
00320
00321 template<typename Pixel>
00322 static Pixel upy(const ImageType & phi, int x, int y)
00323 {
00324 const Pixel u3 = u.neumann_pix2d(x,y+2);
00325 const Pixel u2 = u.neumann_pix2d(x,y+1);
00326 const Pixel u1 = u.neumann_pix2d(x,y);
00327 const Pixel u0 = u.neumann_pix2d(x,y-1);
00328 return combine(u3-u2,u2-u1,u1-u0);
00329 }
00330
00331 template<typename Pixel>
00332 static Pixel upy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00333 {
00334 const Pixel u3 = u.neumann_pix3d(x,y+2,z);
00335 const Pixel u2 = u.neumann_pix3d(x,y+1,z);
00336 const Pixel u1 = u.neumann_pix3d(x,y,z);
00337 const Pixel u0 = u.neumann_pix3d(x,y-1,z);
00338 return combine(u3-u2,u2-u1,u1-u0);
00339 }
00340
00341 template<typename Pixel>
00342 static Pixel umy(const ImageType & phi, int x, int y)
00343 {
00344 const Pixel u3 = u.neumann_pix2d(x,y+1);
00345 const Pixel u2 = u.neumann_pix2d(x,y);
00346 const Pixel u1 = u.neumann_pix2d(x,y-1);
00347 const Pixel u0 = u.neumann_pix2d(x,y-2);
00348 return combine(u1-u0,u2-u1,u3-u2);
00349 }
00350
00351 template<typename Pixel>
00352 static Pixel umy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00353 {
00354 const Pixel u3 = u.neumann_pix3d(x,y+1,z);
00355 const Pixel u2 = u.neumann_pix3d(x,y,z);
00356 const Pixel u1 = u.neumann_pix3d(x,y-1,z);
00357 const Pixel u0 = u.neumann_pix3d(x,y-2,z);
00358 return combine(u1-u0,u2-u1,u3-u2);
00359 }
00360
00361 template<typename Pixel>
00362 static Pixel upz(const MImage<Pixel, 3> &u, int x, int y, int z)
00363 {
00364 const Pixel u3 = u.neumann_pix3d(x,y,z+2);
00365 const Pixel u2 = u.neumann_pix3d(x,y,z+1);
00366 const Pixel u1 = u.neumann_pix3d(x,y,z);
00367 const Pixel u0 = u.neumann_pix3d(x,y,z-1);
00368 return combine(u3-u2,u2-u1,u1-u0);
00369 }
00370
00371 template<typename Pixel>
00372 static Pixel umz(const MImage<Pixel, 3> &u, int x, int y, int z)
00373 {
00374 const Pixel u3 = u.neumann_pix3d(x,y,z+1);
00375 const Pixel u2 = u.neumann_pix3d(x,y,z);
00376 const Pixel u1 = u.neumann_pix3d(x,y,z-1);
00377 const Pixel u0 = u.neumann_pix3d(x,y,z-2);
00378 return combine(u1-u0,u2-u1,u3-u2);
00379 }
00380 };
00381
00382
00387 class Weno5
00388 {
00389 public:
00390
00391 template<typename Pixel, ftk::IndexType VDimension>
00392 static Pixel combine(const Pixel v1, const Pixel v2, const Pixel v3, const Pixel v4, const Pixel v5)
00393 {
00394
00395 Pixel s1 = 13 * sqr(v1 - 2*v2 + v3) / 12 + sqr(v1 - 4*v2 + 3*v3) / 4;
00396 Pixel s2 = 13 * sqr(v2 - 2*v3 + v4) / 12 + sqr(v2 - v4) / 4;
00397 Pixel s3 = 13 * sqr(v3 - 2*v4 + v5) / 12 + sqr(3*v3 - 4*v4 + v5) / 4;
00398
00399 s1 = 1 / sqr(Globals::EPS + s1) / 10;
00400 s2 = 6 / sqr(Globals::EPS + s2) / 10;
00401 s3 = 3 / sqr(Globals::EPS + s3) / 10;
00402
00403 const Pixel s = s1 + s2 + s3;
00404 s1 /= s;
00405 s2 /= s;
00406 s3 /= s;
00407
00408 return ( s1 * (2*v1 - 7*v2 + 11*v3) + s2 * (-v2 + 5*v3 + 2*v4 ) + s3 * (2*v3 + 5*v4 - v5) ) / 6;
00409 }
00410
00411 template<typename Pixel>
00412 static Pixel upx(const ImageType & phi, int x, int y = 0)
00413 {
00414 const Pixel u5 = u.neumann_pix2d(x+3,y);
00415 const Pixel u4 = u.neumann_pix2d(x+2,y);
00416 const Pixel u3 = u.neumann_pix2d(x+1,y);
00417 const Pixel u2 = u.neumann_pix2d(x,y);
00418 const Pixel u1 = u.neumann_pix2d(x-1,y);
00419 const Pixel u0 = u.neumann_pix2d(x-2,y);
00420 return combine(u5-u4,u4-u3,u3-u2,u2-u1,u1-u0);
00421 }
00422
00423 template<typename Pixel>
00424 static Pixel upx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00425 {
00426 const Pixel u5 = u.neumann_pix3d(x+3,y,z);
00427 const Pixel u4 = u.neumann_pix3d(x+2,y,z);
00428 const Pixel u3 = u.neumann_pix3d(x+1,y,z);
00429 const Pixel u2 = u.neumann_pix3d(x,y,z);
00430 const Pixel u1 = u.neumann_pix3d(x-1,y,z);
00431 const Pixel u0 = u.neumann_pix3d(x-2,y,z);
00432 return combine(u5-u4,u4-u3,u3-u2,u2-u1,u1-u0);
00433 }
00434
00435 template<typename Pixel>
00436 static Pixel umx(const ImageType & phi, int x, int y = 0)
00437 {
00438 const Pixel u5 = u.neumann_pix2d(x+2,y);
00439 const Pixel u4 = u.neumann_pix2d(x+1,y);
00440 const Pixel u3 = u.neumann_pix2d(x,y);
00441 const Pixel u2 = u.neumann_pix2d(x-1,y);
00442 const Pixel u1 = u.neumann_pix2d(x-2,y);
00443 const Pixel u0 = u.neumann_pix2d(x-3,y);
00444 return combine(u1-u0,u2-u1,u3-u2,u4-u3,u5-u4);
00445 }
00446
00447 template<typename Pixel>
00448 static Pixel umx(const MImage<Pixel, 3> &u, int x, int y = 0, int z = 0)
00449 {
00450 const Pixel u5 = u.neumann_pix3d(x+2,y,z);
00451 const Pixel u4 = u.neumann_pix3d(x+1,y,z);
00452 const Pixel u3 = u.neumann_pix3d(x,y,z);
00453 const Pixel u2 = u.neumann_pix3d(x-1,y,z);
00454 const Pixel u1 = u.neumann_pix3d(x-2,y,z);
00455 const Pixel u0 = u.neumann_pix3d(x-3,y,z);
00456 return combine(u1-u0,u2-u1,u3-u2,u4-u3,u5-u4);
00457 }
00458
00459 template<typename Pixel>
00460 static Pixel upy(const ImageType & phi, int x, int y, int z = 0)
00461 {
00462 const Pixel u5 = u.neumann_pix2d(x,y+3);
00463 const Pixel u4 = u.neumann_pix2d(x,y+2);
00464 const Pixel u3 = u.neumann_pix2d(x,y+1);
00465 const Pixel u2 = u.neumann_pix2d(x,y);
00466 const Pixel u1 = u.neumann_pix2d(x,y-1);
00467 const Pixel u0 = u.neumann_pix2d(x,y-2);
00468 return combine(u5-u4,u4-u3,u3-u2,u2-u1,u1-u0);
00469 }
00470
00471 template<typename Pixel>
00472 static Pixel upy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00473 {
00474 const Pixel u5 = u.neumann_pix3d(x,y+3,z);
00475 const Pixel u4 = u.neumann_pix3d(x,y+2,z);
00476 const Pixel u3 = u.neumann_pix3d(x,y+1,z);
00477 const Pixel u2 = u.neumann_pix3d(x,y,z);
00478 const Pixel u1 = u.neumann_pix3d(x,y-1,z);
00479 const Pixel u0 = u.neumann_pix3d(x,y-2,z);
00480 return combine(u5-u4,u4-u3,u3-u2,u2-u1,u1-u0);
00481 }
00482
00483 template<typename Pixel>
00484 static Pixel umy(const ImageType & phi, int x, int y, int z = 0)
00485 {
00486 const Pixel u5 = u.neumann_pix2d(x,y+2);
00487 const Pixel u4 = u.neumann_pix2d(x,y+1);
00488 const Pixel u3 = u.neumann_pix2d(x,y);
00489 const Pixel u2 = u.neumann_pix2d(x,y-1);
00490 const Pixel u1 = u.neumann_pix2d(x,y-2);
00491 const Pixel u0 = u.neumann_pix2d(x,y-3);
00492 return combine(u1-u0,u2-u1,u3-u2,u4-u3,u5-u4);
00493 }
00494
00495 template<typename Pixel>
00496 static Pixel umy(const MImage<Pixel, 3> &u, int x, int y, int z = 0)
00497 {
00498 const Pixel u5 = u.neumann_pix3d(x,y+2,z);
00499 const Pixel u4 = u.neumann_pix3d(x,y+1,z);
00500 const Pixel u3 = u.neumann_pix3d(x,y,z);
00501 const Pixel u2 = u.neumann_pix3d(x,y-1,z);
00502 const Pixel u1 = u.neumann_pix3d(x,y-2,z);
00503 const Pixel u0 = u.neumann_pix3d(x,y-3,z);
00504 return combine(u1-u0,u2-u1,u3-u2,u4-u3,u5-u4);
00505 }
00506
00507 template<typename Pixel>
00508 static Pixel upz(const MImage<Pixel, 3> &u, int x, int y, int z)
00509 {
00510 const Pixel u5 = u.neumann_pix3d(x,y,z+3);
00511 const Pixel u4 = u.neumann_pix3d(x,y,z+2);
00512 const Pixel u3 = u.neumann_pix3d(x,y,z+1);
00513 const Pixel u2 = u.neumann_pix3d(x,y,z);
00514 const Pixel u1 = u.neumann_pix3d(x,y,z-1);
00515 const Pixel u0 = u.neumann_pix3d(x,y,z-2);
00516 return combine(u5-u4,u4-u3,u3-u2,u2-u1,u1-u0);
00517 }
00518
00519 template<typename Pixel>
00520 static Pixel umz(const MImage<Pixel, 3> &u, int x, int y, int z)
00521 {
00522 const Pixel u5 = u.neumann_pix3d(x,y,z+2);
00523 const Pixel u4 = u.neumann_pix3d(x,y,z+1);
00524 const Pixel u3 = u.neumann_pix3d(x,y,z);
00525 const Pixel u2 = u.neumann_pix3d(x,y,z-1);
00526 const Pixel u1 = u.neumann_pix3d(x,y,z-2);
00527 const Pixel u0 = u.neumann_pix3d(x,y,z-3);
00528 return combine(u1-u0,u2-u1,u3-u2,u4-u3,u5-u4);
00529 }
00530 };
00531
00534 #endif
00535
00536 }
00537
00538 #endif //APPROXIMATORS_H