45 #ifndef CLIPPER_MAP_INTERP
46 #define CLIPPER_MAP_INTERP
70 template<
class T,
class M>
static void interp(
const M& map,
const Coord_map& pos, T& val );
71 inline static int order() {
return 0; }
91 template<
class T,
class M>
static void interp(
const M& map,
const Coord_map& pos, T& val );
92 inline static int order() {
return 1; }
112 template<
class T,
class M>
static void interp(
const M& map,
const Coord_map& pos, T& val );
115 inline static int order() {
return 3; }
152 if ( map.in_map( c ) ) {
154 return map.in_map( c );
168 ftype u0 = floor( pos.
u() );
169 ftype v0 = floor( pos.
v() );
170 ftype w0 = floor( pos.
w() );
171 typename M::Map_reference_coord
172 r( map,
Coord_grid(
int(u0),
int(v0),
int(w0) ) );
173 T cu1( pos.
u() - u0 );
174 T cv1( pos.
v() - v0 );
175 T cw1( pos.
w() - w0 );
179 T r00 = cw0 * map[ r ];
180 r00 += cw1 * map[ r.next_w() ];
181 T r01 = cw1 * map[ r.next_v() ];
182 r01 += cw0 * map[ r.prev_w() ];
183 T r11 = cw0 * map[ r.next_u() ];
184 r11 += cw1 * map[ r.next_w() ];
185 T r10 = cw1 * map[ r.prev_v() ];
186 r10 += cw0 * map[ r.prev_w() ];
187 val = ( cu0*( cv0*r00 + cv1*r01 ) + cu1*( cv0*r10 + cv1*r11 ) );
201 if ( map.in_map( c ) ) {
203 return map.in_map( c );
215 ftype u0 = floor( pos.
u() );
216 ftype v0 = floor( pos.
v() );
217 ftype w0 = floor( pos.
w() );
218 typename M::Map_reference_coord iw, iv,
219 iu( map,
Coord_grid(
int(u0)-1,
int(v0)-1,
int(w0)-1 ) );
220 T su, sv, sw, cu[4], cv[4], cw[4];
221 T cu1( pos.
u() - u0 );
222 T cv1( pos.
v() - v0 );
223 T cw1( pos.
w() - w0 );
227 cu[0] = -0.5*cu1*cu0*cu0;
228 cu[1] = cu0*( -1.5*cu1*cu1 + cu1 + 1.0 );
229 cu[2] = cu1*( -1.5*cu0*cu0 + cu0 + 1.0 );
230 cu[3] = -0.5*cu1*cu1*cu0;
231 cv[0] = -0.5*cv1*cv0*cv0;
232 cv[1] = cv0*( -1.5*cv1*cv1 + cv1 + 1.0 );
233 cv[2] = cv1*( -1.5*cv0*cv0 + cv0 + 1.0 );
234 cv[3] = -0.5*cv1*cv1*cv0;
235 cw[0] = -0.5*cw1*cw0*cw0;
236 cw[1] = cw0*( -1.5*cw1*cw1 + cw1 + 1.0 );
237 cw[2] = cw1*( -1.5*cw0*cw0 + cw0 + 1.0 );
238 cw[3] = -0.5*cw1*cw1*cw0;
241 for ( j = 0; j < 4; j++ ) {
244 for ( i = 0; i < 4; i++ ) {
246 sw = cw[0] * T( map[ iw ] );
247 sw += cw[1] * T( map[ iw.next_w() ] );
248 sw += cw[2] * T( map[ iw.next_w() ] );
249 sw += cw[3] * T( map[ iw.next_w() ] );
269 ftype u0 = floor( pos.
u() );
270 ftype v0 = floor( pos.
v() );
271 ftype w0 = floor( pos.
w() );
272 typename M::Map_reference_coord iw, iv,
273 iu( map,
Coord_grid(
int(u0)-1,
int(v0)-1,
int(w0)-1 ) );
274 T s1, s2, s3, du1, dv1, dv2, dw1, dw2, dw3;
275 T cu[4], cv[4], cw[4], gu[4], gv[4], gw[4];
276 T cu1( pos.
u() - u0 );
277 T cv1( pos.
v() - v0 );
278 T cw1( pos.
w() - w0 );
282 cu[0] = -0.5*cu1*cu0*cu0;
283 cu[1] = cu0*( -1.5*cu1*cu1 + cu1 + 1.0 );
284 cu[2] = cu1*( -1.5*cu0*cu0 + cu0 + 1.0 );
285 cu[3] = -0.5*cu1*cu1*cu0;
286 cv[0] = -0.5*cv1*cv0*cv0;
287 cv[1] = cv0*( -1.5*cv1*cv1 + cv1 + 1.0 );
288 cv[2] = cv1*( -1.5*cv0*cv0 + cv0 + 1.0 );
289 cv[3] = -0.5*cv1*cv1*cv0;
290 cw[0] = -0.5*cw1*cw0*cw0;
291 cw[1] = cw0*( -1.5*cw1*cw1 + cw1 + 1.0 );
292 cw[2] = cw1*( -1.5*cw0*cw0 + cw0 + 1.0 );
293 cw[3] = -0.5*cw1*cw1*cw0;
294 gu[0] = cu0*( 1.5*cu1 - 0.5 );
295 gu[1] = cu1*( 4.5*cu1 - 5.0 );
296 gu[2] = -cu0*( 4.5*cu0 - 5.0 );
297 gu[3] = -cu1*( 1.5*cu0 - 0.5 );
298 gv[0] = cv0*( 1.5*cv1 - 0.5 );
299 gv[1] = cv1*( 4.5*cv1 - 5.0 );
300 gv[2] = -cv0*( 4.5*cv0 - 5.0 );
301 gv[3] = -cv1*( 1.5*cv0 - 0.5 );
302 gw[0] = cw0*( 1.5*cw1 - 0.5 );
303 gw[1] = cw1*( 4.5*cw1 - 5.0 );
304 gw[2] = -cw0*( 4.5*cw0 - 5.0 );
305 gw[3] = -cw1*( 1.5*cw0 - 0.5 );
306 s1 = du1 = dv1 = dw1 = 0.0;
308 for ( j = 0; j < 4; j++ ) {
310 s2 = dv2 = dw2 = 0.0;
311 for ( i = 0; i < 4; i++ ) {
313 s3 = cw[0] * T( map[ iw ] );
314 dw3 = gw[0] * T( map[ iw ] );
316 s3 += cw[1] * T( map[ iw ] );
317 dw3 += gw[1] * T( map[ iw ] );
319 s3 += cw[2] * T( map[ iw ] );
320 dw3 += gw[2] * T( map[ iw ] );
322 s3 += cw[3] * T( map[ iw ] );
323 dw3 += gw[3] * T( map[ iw ] );
349 ftype u0 = floor( pos.
u() );
350 ftype v0 = floor( pos.
v() );
351 ftype w0 = floor( pos.
w() );
352 typename M::Map_reference_coord iw, iv,
353 iu( map,
Coord_grid(
int(u0)-1,
int(v0)-1,
int(w0)-1 ) );
354 T s1, s2, s3, du1, dv1, dv2, dw1, dw2, dw3;
355 T duv1, duw1, dvw1, dvw2, duu1, dvv1, dvv2, dww1, dww2, dww3;
356 T cu[4], cv[4], cw[4], gu[4], gv[4], gw[4], ggu[4], ggv[4], ggw[4];
357 T cu1( pos.
u() - u0 );
358 T cv1( pos.
v() - v0 );
359 T cw1( pos.
w() - w0 );
363 cu[0] = -0.5*cu1*cu0*cu0;
364 cu[1] = cu0*( -1.5*cu1*cu1 + cu1 + 1.0 );
365 cu[2] = cu1*( -1.5*cu0*cu0 + cu0 + 1.0 );
366 cu[3] = -0.5*cu1*cu1*cu0;
367 cv[0] = -0.5*cv1*cv0*cv0;
368 cv[1] = cv0*( -1.5*cv1*cv1 + cv1 + 1.0 );
369 cv[2] = cv1*( -1.5*cv0*cv0 + cv0 + 1.0 );
370 cv[3] = -0.5*cv1*cv1*cv0;
371 cw[0] = -0.5*cw1*cw0*cw0;
372 cw[1] = cw0*( -1.5*cw1*cw1 + cw1 + 1.0 );
373 cw[2] = cw1*( -1.5*cw0*cw0 + cw0 + 1.0 );
374 cw[3] = -0.5*cw1*cw1*cw0;
375 gu[0] = cu0*( 1.5*cu1 - 0.5 );
376 gu[1] = cu1*( 4.5*cu1 - 5.0 );
377 gu[2] = -cu0*( 4.5*cu0 - 5.0 );
378 gu[3] = -cu1*( 1.5*cu0 - 0.5 );
379 gv[0] = cv0*( 1.5*cv1 - 0.5 );
380 gv[1] = cv1*( 4.5*cv1 - 5.0 );
381 gv[2] = -cv0*( 4.5*cv0 - 5.0 );
382 gv[3] = -cv1*( 1.5*cv0 - 0.5 );
383 gw[0] = cw0*( 1.5*cw1 - 0.5 );
384 gw[1] = cw1*( 4.5*cw1 - 5.0 );
385 gw[2] = -cw0*( 4.5*cw0 - 5.0 );
386 gw[3] = -cw1*( 1.5*cw0 - 0.5 );
387 ggu[0] = 2.0 - 3.0*cu1;
388 ggu[1] = 9.0*cu1 - 5.0;
389 ggu[2] = 9.0*cu0 - 5.0;
390 ggu[3] = 2.0 - 3.0*cu0;
391 ggv[0] = 2.0 - 3.0*cv1;
392 ggv[1] = 9.0*cv1 - 5.0;
393 ggv[2] = 9.0*cv0 - 5.0;
394 ggv[3] = 2.0 - 3.0*cv0;
395 ggw[0] = 2.0 - 3.0*cw1;
396 ggw[1] = 9.0*cw1 - 5.0;
397 ggw[2] = 9.0*cw0 - 5.0;
398 ggw[3] = 2.0 - 3.0*cw0;
399 s1 = du1 = dv1 = dw1 = duv1 = duw1 = dvw1 = duu1 = dvv1 = dww1 = 0.0;
401 for ( j = 0; j < 4; j++ ) {
403 s2 = dv2 = dw2 = dvw2 = dvv2 = dww2 = 0.0;
404 for ( i = 0; i < 4; i++ ) {
406 s3 = cw[0] * T( map[ iw ] );
407 dw3 = gw[0] * T( map[ iw ] );
408 dww3 = ggw[0] * T( map[ iw ] );
410 s3 += cw[1] * T( map[ iw ] );
411 dw3 += gw[1] * T( map[ iw ] );
412 dww3 += ggw[1] * T( map[ iw ] );
414 s3 += cw[2] * T( map[ iw ] );
415 dw3 += gw[2] * T( map[ iw ] );
416 dww3 += ggw[2] * T( map[ iw ] );
418 s3 += cw[3] * T( map[ iw ] );
419 dw3 += gw[3] * T( map[ iw ] );
420 dww3 += ggw[3] * T( map[ iw ] );
426 dww2 += cv[i] * dww3;
435 dvw1 += cu[j] * dvw2;
437 dvv1 += cu[j] * dvv2;
438 dww1 += cu[j] * dww2;
445 duw1, dvw1, dww1 ) );
const int & u() const
get u
Definition: coords.h:248
static int order()
Order of interpolant.
Definition: map_interp.h:71
static void interp(const M &map, const Coord_map &pos, T &val)
Interpolate map M using type T at coord.
Definition: map_interp.h:166
static void interp(const M &map, const Coord_map &pos, T &val)
Interpolate map M using type T at coord.
Definition: map_interp.h:213
static bool can_interp(const M &map, const Coord_map &pos)
Test if we can interpolate in map M at coord.
Definition: map_interp.h:197
map coordinate: this is like Coord_grid, but non-integer
Definition: coords.h:387
Grid coordinate.
Definition: coords.h:236
Coord_grid floor() const
return integer Coord_grid below this coordinate
Definition: coords.h:405
const ftype & u() const
get u
Definition: coords.h:408
static void interp(const M &map, const Coord_map &pos, T &val)
Interpolate map M using type T at coord.
Definition: map_interp.h:138
Wrapper class for zeroth-order (nearest neighbour) interpolation fns.
Definition: map_interp.h:66
static bool can_interp(const M &map, const Coord_map &pos)
Test if we can interpolate in map M at coord.
Definition: map_interp.h:148
const ftype & w() const
get w
Definition: coords.h:410
ftype64 ftype
ftype definition for floating point representation
Definition: clipper_precision.h:58
Wrapper class for third-order (cubic) interpolation fns.
Definition: map_interp.h:108
static void interp_curv(const M &map, const Coord_map &pos, T &val, Grad_map< T > &grad, Curv_map< T > &curv)
Definition: map_interp.h:347
map coordinate curvatures, with respect to grid u,v,w
Definition: derivs.h:59
Wrapper class for first-order (linear) interpolation fns.
Definition: map_interp.h:87
static int order()
Order of interpolant.
Definition: map_interp.h:92
static void interp_grad(const M &map, const Coord_map &pos, T &val, Grad_map< T > &grad)
Definition: map_interp.h:267
Coord_grid coord_grid() const
return integer Coord_grid nearest this coordinate
Definition: coords.h:403
map coordinate gradient, with respect to grid u,v,w
Definition: derivs.h:56
const ftype & v() const
get v
Definition: coords.h:409
static bool can_interp(const M &map, const Coord_map &pos)
Test if we can interpolate in map M at coord.
Definition: map_interp.h:128
static int order()
Order of interpolant.
Definition: map_interp.h:115
3x3-matrix class
Definition: clipper_types.h:182