00001
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #ifndef CLIPPER_MAP_INTERP
00046 #define CLIPPER_MAP_INTERP
00047
00048 #include "derivs.h"
00049
00050
00051 namespace clipper
00052 {
00053
00055
00066 class Interp_nearest
00067 {
00068 public:
00069 template<class M> static bool can_interp( const M& map, const Coord_map& pos );
00070 template<class T, class M> static void interp( const M& map, const Coord_map& pos, T& val );
00071 inline static int order() { return 0; }
00072 };
00073
00075
00087 class Interp_linear
00088 {
00089 public:
00090 template<class M> static bool can_interp( const M& map, const Coord_map& pos );
00091 template<class T, class M> static void interp( const M& map, const Coord_map& pos, T& val );
00092 inline static int order() { return 1; }
00093 };
00094
00096
00108 class Interp_cubic
00109 {
00110 public:
00111 template<class M> static bool can_interp( const M& map, const Coord_map& pos );
00112 template<class T, class M> static void interp( const M& map, const Coord_map& pos, T& val );
00113 template<class T, class M> static void interp_grad( const M& map, const Coord_map& pos, T& val, Grad_map<T>& grad );
00114 template<class T, class M> static void interp_curv( const M& map, const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv );
00115 inline static int order() { return 3; }
00116 };
00117
00118
00119
00120
00121
00128 template<class M> bool Interp_nearest::can_interp( const M& map, const Coord_map& pos )
00129 { return map.in_map( pos.coord_grid() ); }
00130
00138 template<class T, class M> void Interp_nearest::interp( const M& map, const Coord_map& pos, T& val )
00139 { val = map.get_data( pos.coord_grid() ); }
00140
00141
00148 template<class M> bool Interp_linear::can_interp( const M& map, const Coord_map& pos )
00149 {
00150 Coord_grid c( pos.floor() );
00151 c.u() -= order()/2; c.v() -= order()/2; c.w() -= order()/2;
00152 if ( map.in_map( c ) ) {
00153 c.u() += order(); c.v() += order(); c.w() += order();
00154 return map.in_map( c );
00155 }
00156 return false;
00157 }
00158
00166 template<class T, class M> void Interp_linear::interp( const M& map, const Coord_map& pos, T& val )
00167 {
00168 ftype u0 = floor( pos.u() );
00169 ftype v0 = floor( pos.v() );
00170 ftype w0 = floor( pos.w() );
00171 typename M::Map_reference_coord
00172 r( map, Coord_grid( int(u0), int(v0), int(w0) ) );
00173 T cu1( pos.u() - u0 );
00174 T cv1( pos.v() - v0 );
00175 T cw1( pos.w() - w0 );
00176 T cu0( 1.0 - cu1 );
00177 T cv0( 1.0 - cv1 );
00178 T cw0( 1.0 - cw1 );
00179 T r00 = cw0 * map[ r ];
00180 r00 += cw1 * map[ r.next_w() ];
00181 T r01 = cw1 * map[ r.next_v() ];
00182 r01 += cw0 * map[ r.prev_w() ];
00183 T r11 = cw0 * map[ r.next_u() ];
00184 r11 += cw1 * map[ r.next_w() ];
00185 T r10 = cw1 * map[ r.prev_v() ];
00186 r10 += cw0 * map[ r.prev_w() ];
00187 val = ( cu0*( cv0*r00 + cv1*r01 ) + cu1*( cv0*r10 + cv1*r11 ) );
00188 }
00189
00190
00197 template<class M> bool Interp_cubic::can_interp( const M& map, const Coord_map& pos )
00198 {
00199 Coord_grid c( pos.floor() );
00200 c.u() -= order()/2; c.v() -= order()/2; c.w() -= order()/2;
00201 if ( map.in_map( c ) ) {
00202 c.u() += order(); c.v() += order(); c.w() += order();
00203 return map.in_map( c );
00204 }
00205 return false;
00206 }
00207
00213 template<class T, class M> void Interp_cubic::interp( const M& map, const Coord_map& pos, T& val )
00214 {
00215 ftype u0 = floor( pos.u() );
00216 ftype v0 = floor( pos.v() );
00217 ftype w0 = floor( pos.w() );
00218 typename M::Map_reference_coord iw, iv,
00219 iu( map, Coord_grid( int(u0)-1, int(v0)-1, int(w0)-1 ) );
00220 T su, sv, sw, cu[4], cv[4], cw[4];
00221 T cu1( pos.u() - u0 );
00222 T cv1( pos.v() - v0 );
00223 T cw1( pos.w() - w0 );
00224 T cu0( 1.0 - cu1 );
00225 T cv0( 1.0 - cv1 );
00226 T cw0( 1.0 - cw1 );
00227 cu[0] = -0.5*cu1*cu0*cu0;
00228 cu[1] = cu0*( -1.5*cu1*cu1 + cu1 + 1.0 );
00229 cu[2] = cu1*( -1.5*cu0*cu0 + cu0 + 1.0 );
00230 cu[3] = -0.5*cu1*cu1*cu0;
00231 cv[0] = -0.5*cv1*cv0*cv0;
00232 cv[1] = cv0*( -1.5*cv1*cv1 + cv1 + 1.0 );
00233 cv[2] = cv1*( -1.5*cv0*cv0 + cv0 + 1.0 );
00234 cv[3] = -0.5*cv1*cv1*cv0;
00235 cw[0] = -0.5*cw1*cw0*cw0;
00236 cw[1] = cw0*( -1.5*cw1*cw1 + cw1 + 1.0 );
00237 cw[2] = cw1*( -1.5*cw0*cw0 + cw0 + 1.0 );
00238 cw[3] = -0.5*cw1*cw1*cw0;
00239 su = 0.0;
00240 int i, j;
00241 for ( j = 0; j < 4; j++ ) {
00242 iv = iu;
00243 sv = 0.0;
00244 for ( i = 0; i < 4; i++ ) {
00245 iw = iv;
00246 sw = cw[0] * T( map[ iw ] );
00247 sw += cw[1] * T( map[ iw.next_w() ] );
00248 sw += cw[2] * T( map[ iw.next_w() ] );
00249 sw += cw[3] * T( map[ iw.next_w() ] );
00250 sv += cv[i] * sw;
00251 iv.next_v();
00252 }
00253 su += cu[j] * sv;
00254 iu.next_u();
00255 }
00256 val = su;
00257 }
00258
00259
00267 template<class T, class M> void Interp_cubic::interp_grad( const M& map, const Coord_map& pos, T& val, Grad_map<T>& grad )
00268 {
00269 ftype u0 = floor( pos.u() );
00270 ftype v0 = floor( pos.v() );
00271 ftype w0 = floor( pos.w() );
00272 typename M::Map_reference_coord iw, iv,
00273 iu( map, Coord_grid( int(u0)-1, int(v0)-1, int(w0)-1 ) );
00274 T s1, s2, s3, du1, dv1, dv2, dw1, dw2, dw3;
00275 T cu[4], cv[4], cw[4], gu[4], gv[4], gw[4];
00276 T cu1( pos.u() - u0 );
00277 T cv1( pos.v() - v0 );
00278 T cw1( pos.w() - w0 );
00279 T cu0( 1.0 - cu1 );
00280 T cv0( 1.0 - cv1 );
00281 T cw0( 1.0 - cw1 );
00282 cu[0] = -0.5*cu1*cu0*cu0;
00283 cu[1] = cu0*( -1.5*cu1*cu1 + cu1 + 1.0 );
00284 cu[2] = cu1*( -1.5*cu0*cu0 + cu0 + 1.0 );
00285 cu[3] = -0.5*cu1*cu1*cu0;
00286 cv[0] = -0.5*cv1*cv0*cv0;
00287 cv[1] = cv0*( -1.5*cv1*cv1 + cv1 + 1.0 );
00288 cv[2] = cv1*( -1.5*cv0*cv0 + cv0 + 1.0 );
00289 cv[3] = -0.5*cv1*cv1*cv0;
00290 cw[0] = -0.5*cw1*cw0*cw0;
00291 cw[1] = cw0*( -1.5*cw1*cw1 + cw1 + 1.0 );
00292 cw[2] = cw1*( -1.5*cw0*cw0 + cw0 + 1.0 );
00293 cw[3] = -0.5*cw1*cw1*cw0;
00294 gu[0] = cu0*( 1.5*cu1 - 0.5 );
00295 gu[1] = cu1*( 4.5*cu1 - 5.0 );
00296 gu[2] = -cu0*( 4.5*cu0 - 5.0 );
00297 gu[3] = -cu1*( 1.5*cu0 - 0.5 );
00298 gv[0] = cv0*( 1.5*cv1 - 0.5 );
00299 gv[1] = cv1*( 4.5*cv1 - 5.0 );
00300 gv[2] = -cv0*( 4.5*cv0 - 5.0 );
00301 gv[3] = -cv1*( 1.5*cv0 - 0.5 );
00302 gw[0] = cw0*( 1.5*cw1 - 0.5 );
00303 gw[1] = cw1*( 4.5*cw1 - 5.0 );
00304 gw[2] = -cw0*( 4.5*cw0 - 5.0 );
00305 gw[3] = -cw1*( 1.5*cw0 - 0.5 );
00306 s1 = du1 = dv1 = dw1 = 0.0;
00307 int i, j;
00308 for ( j = 0; j < 4; j++ ) {
00309 iv = iu;
00310 s2 = dv2 = dw2 = 0.0;
00311 for ( i = 0; i < 4; i++ ) {
00312 iw = iv;
00313 s3 = cw[0] * T( map[ iw ] );
00314 dw3 = gw[0] * T( map[ iw ] );
00315 iw.next_w();
00316 s3 += cw[1] * T( map[ iw ] );
00317 dw3 += gw[1] * T( map[ iw ] );
00318 iw.next_w();
00319 s3 += cw[2] * T( map[ iw ] );
00320 dw3 += gw[2] * T( map[ iw ] );
00321 iw.next_w();
00322 s3 += cw[3] * T( map[ iw ] );
00323 dw3 += gw[3] * T( map[ iw ] );
00324 s2 += cv[i] * s3;
00325 dv2 += gv[i] * s3;
00326 dw2 += cv[i] * dw3;
00327 iv.next_v();
00328 }
00329 s1 += cu[j] * s2;
00330 du1 += gu[j] * s2;
00331 dv1 += cu[j] * dv2;
00332 dw1 += cu[j] * dw2;
00333 iu.next_u();
00334 }
00335 val = s1;
00336 grad = Grad_map<T>( du1, dv1, dw1 );
00337 }
00338
00339
00347 template<class T, class M> void Interp_cubic::interp_curv( const M& map, const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv )
00348 {
00349 ftype u0 = floor( pos.u() );
00350 ftype v0 = floor( pos.v() );
00351 ftype w0 = floor( pos.w() );
00352 typename M::Map_reference_coord iw, iv,
00353 iu( map, Coord_grid( int(u0)-1, int(v0)-1, int(w0)-1 ) );
00354 T s1, s2, s3, du1, dv1, dv2, dw1, dw2, dw3;
00355 T duv1, duw1, dvw1, dvw2, duu1, dvv1, dvv2, dww1, dww2, dww3;
00356 T cu[4], cv[4], cw[4], gu[4], gv[4], gw[4], ggu[4], ggv[4], ggw[4];
00357 T cu1( pos.u() - u0 );
00358 T cv1( pos.v() - v0 );
00359 T cw1( pos.w() - w0 );
00360 T cu0( 1.0 - cu1 );
00361 T cv0( 1.0 - cv1 );
00362 T cw0( 1.0 - cw1 );
00363 cu[0] = -0.5*cu1*cu0*cu0;
00364 cu[1] = cu0*( -1.5*cu1*cu1 + cu1 + 1.0 );
00365 cu[2] = cu1*( -1.5*cu0*cu0 + cu0 + 1.0 );
00366 cu[3] = -0.5*cu1*cu1*cu0;
00367 cv[0] = -0.5*cv1*cv0*cv0;
00368 cv[1] = cv0*( -1.5*cv1*cv1 + cv1 + 1.0 );
00369 cv[2] = cv1*( -1.5*cv0*cv0 + cv0 + 1.0 );
00370 cv[3] = -0.5*cv1*cv1*cv0;
00371 cw[0] = -0.5*cw1*cw0*cw0;
00372 cw[1] = cw0*( -1.5*cw1*cw1 + cw1 + 1.0 );
00373 cw[2] = cw1*( -1.5*cw0*cw0 + cw0 + 1.0 );
00374 cw[3] = -0.5*cw1*cw1*cw0;
00375 gu[0] = cu0*( 1.5*cu1 - 0.5 );
00376 gu[1] = cu1*( 4.5*cu1 - 5.0 );
00377 gu[2] = -cu0*( 4.5*cu0 - 5.0 );
00378 gu[3] = -cu1*( 1.5*cu0 - 0.5 );
00379 gv[0] = cv0*( 1.5*cv1 - 0.5 );
00380 gv[1] = cv1*( 4.5*cv1 - 5.0 );
00381 gv[2] = -cv0*( 4.5*cv0 - 5.0 );
00382 gv[3] = -cv1*( 1.5*cv0 - 0.5 );
00383 gw[0] = cw0*( 1.5*cw1 - 0.5 );
00384 gw[1] = cw1*( 4.5*cw1 - 5.0 );
00385 gw[2] = -cw0*( 4.5*cw0 - 5.0 );
00386 gw[3] = -cw1*( 1.5*cw0 - 0.5 );
00387 ggu[0] = 2.0 - 3.0*cu1;
00388 ggu[1] = 9.0*cu1 - 5.0;
00389 ggu[2] = 9.0*cu0 - 5.0;
00390 ggu[3] = 2.0 - 3.0*cu0;
00391 ggv[0] = 2.0 - 3.0*cv1;
00392 ggv[1] = 9.0*cv1 - 5.0;
00393 ggv[2] = 9.0*cv0 - 5.0;
00394 ggv[3] = 2.0 - 3.0*cv0;
00395 ggw[0] = 2.0 - 3.0*cw1;
00396 ggw[1] = 9.0*cw1 - 5.0;
00397 ggw[2] = 9.0*cw0 - 5.0;
00398 ggw[3] = 2.0 - 3.0*cw0;
00399 s1 = du1 = dv1 = dw1 = duv1 = duw1 = dvw1 = duu1 = dvv1 = dww1 = 0.0;
00400 int i, j;
00401 for ( j = 0; j < 4; j++ ) {
00402 iv = iu;
00403 s2 = dv2 = dw2 = dvw2 = dvv2 = dww2 = 0.0;
00404 for ( i = 0; i < 4; i++ ) {
00405 iw = iv;
00406 s3 = cw[0] * T( map[ iw ] );
00407 dw3 = gw[0] * T( map[ iw ] );
00408 dww3 = ggw[0] * T( map[ iw ] );
00409 iw.next_w();
00410 s3 += cw[1] * T( map[ iw ] );
00411 dw3 += gw[1] * T( map[ iw ] );
00412 dww3 += ggw[1] * T( map[ iw ] );
00413 iw.next_w();
00414 s3 += cw[2] * T( map[ iw ] );
00415 dw3 += gw[2] * T( map[ iw ] );
00416 dww3 += ggw[2] * T( map[ iw ] );
00417 iw.next_w();
00418 s3 += cw[3] * T( map[ iw ] );
00419 dw3 += gw[3] * T( map[ iw ] );
00420 dww3 += ggw[3] * T( map[ iw ] );
00421 s2 += cv[i] * s3;
00422 dv2 += gv[i] * s3;
00423 dw2 += cv[i] * dw3;
00424 dvw2 += gv[i] * dw3;
00425 dvv2 += ggv[i] * s3;
00426 dww2 += cv[i] * dww3;
00427 iv.next_v();
00428 }
00429 s1 += cu[j] * s2;
00430 du1 += gu[j] * s2;
00431 dv1 += cu[j] * dv2;
00432 dw1 += cu[j] * dw2;
00433 duv1 += gu[j] * dv2;
00434 duw1 += gu[j] * dw2;
00435 dvw1 += cu[j] * dvw2;
00436 duu1 += ggu[j] * s2;
00437 dvv1 += cu[j] * dvv2;
00438 dww1 += cu[j] * dww2;
00439 iu.next_u();
00440 }
00441 val = s1;
00442 grad = Grad_map<T>( du1, dv1, dw1 );
00443 curv = Curv_map<T>( Mat33<T>( duu1, duv1, duw1,
00444 duv1, dvv1, dvw1,
00445 duw1, dvw1, dww1 ) );
00446 }
00447
00448
00449 }
00450
00451
00452 #endif