Clipper
map_interp.h
1
4//C Copyright (C) 2000-2006 Kevin Cowtan and University of York
5//L
6//L This library is free software and is distributed under the terms
7//L and conditions of version 2.1 of the GNU Lesser General Public
8//L Licence (LGPL) with the following additional clause:
9//L
10//L `You may also combine or link a "work that uses the Library" to
11//L produce a work containing portions of the Library, and distribute
12//L that work under terms of your choice, provided that you give
13//L prominent notice with each copy of the work that the specified
14//L version of the Library is used in it, and that you include or
15//L provide public access to the complete corresponding
16//L machine-readable source code for the Library including whatever
17//L changes were used in the work. (i.e. If you make changes to the
18//L Library you must distribute those, but you do not need to
19//L distribute source or object code to those portions of the work
20//L not covered by this licence.)'
21//L
22//L Note that this clause grants an additional right and does not impose
23//L any additional restriction, and so does not affect compatibility
24//L with the GNU General Public Licence (GPL). If you wish to negotiate
25//L other terms, please contact the maintainer.
26//L
27//L You can redistribute it and/or modify the library under the terms of
28//L the GNU Lesser General Public License as published by the Free Software
29//L Foundation; either version 2.1 of the License, or (at your option) any
30//L later version.
31//L
32//L This library is distributed in the hope that it will be useful, but
33//L WITHOUT ANY WARRANTY; without even the implied warranty of
34//L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35//L Lesser General Public License for more details.
36//L
37//L You should have received a copy of the CCP4 licence and/or GNU
38//L Lesser General Public License along with this library; if not, write
39//L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
40//L The GNU Lesser General Public can also be obtained by writing to the
41//L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
42//L MA 02111-1307 USA
43
44
45#ifndef CLIPPER_MAP_INTERP
46#define CLIPPER_MAP_INTERP
47
48#include "derivs.h"
49
50
51namespace clipper
52{
53
55
67 {
68 public:
69 template<class M> static bool can_interp( const M& map, const Coord_map& pos );
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; }
72 };
73
75
88 {
89 public:
90 template<class M> static bool can_interp( const M& map, const Coord_map& pos );
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; }
93 };
94
96
109 {
110 public:
111 template<class M> static bool can_interp( const M& map, const Coord_map& pos );
112 template<class T, class M> static void interp( const M& map, const Coord_map& pos, T& val );
113 template<class T, class M> static void interp_grad( const M& map, const Coord_map& pos, T& val, Grad_map<T>& grad );
114 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 );
115 inline static int order() { return 3; }
116 };
117
118
119
120 // template implementations
121
128 template<class M> bool Interp_nearest::can_interp( const M& map, const Coord_map& pos )
129 { return map.in_map( pos.coord_grid() ); }
130
138 template<class T, class M> void Interp_nearest::interp( const M& map, const Coord_map& pos, T& val )
139 { val = map.get_data( pos.coord_grid() ); }
140
141
148 template<class M> bool Interp_linear::can_interp( const M& map, const Coord_map& pos )
149 {
150 Coord_grid c( pos.floor() ); // for even order change floor to coord_grid
151 c.u() -= order()/2; c.v() -= order()/2; c.w() -= order()/2;
152 if ( map.in_map( c ) ) {
153 c.u() += order(); c.v() += order(); c.w() += order();
154 return map.in_map( c );
155 }
156 return false;
157 }
158
166 template<class T, class M> void Interp_linear::interp( const M& map, const Coord_map& pos, T& val )
167 {
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 );
176 T cu0( 1.0 - cu1 );
177 T cv0( 1.0 - cv1 );
178 T cw0( 1.0 - cw1 );
179 T r00 = cw0 * map[ r ]; // careful with evaluation order
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 ) );
188 }
189
190
197 template<class M> bool Interp_cubic::can_interp( const M& map, const Coord_map& pos )
198 {
199 Coord_grid c( pos.floor() ); // for even order change floor to coord_grid
200 c.u() -= order()/2; c.v() -= order()/2; c.w() -= order()/2;
201 if ( map.in_map( c ) ) {
202 c.u() += order(); c.v() += order(); c.w() += order();
203 return map.in_map( c );
204 }
205 return false;
206 }
207
213 template<class T, class M> void Interp_cubic::interp( const M& map, const Coord_map& pos, T& val )
214 {
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 );
224 T cu0( 1.0 - cu1 );
225 T cv0( 1.0 - cv1 );
226 T cw0( 1.0 - cw1 );
227 cu[0] = -0.5*cu1*cu0*cu0; // cubic spline coeffs: u
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; // cubic spline coeffs: v
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; // cubic spline coeffs: w
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;
239 su = 0.0;
240 int i, j;
241 for ( j = 0; j < 4; j++ ) {
242 iv = iu;
243 sv = 0.0;
244 for ( i = 0; i < 4; i++ ) {
245 iw = iv;
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() ] );
250 sv += cv[i] * sw;
251 iv.next_v();
252 }
253 su += cu[j] * sv;
254 iu.next_u();
255 }
256 val = su;
257 }
258
259
267 template<class T, class M> void Interp_cubic::interp_grad( const M& map, const Coord_map& pos, T& val, Grad_map<T>& grad )
268 {
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 );
279 T cu0( 1.0 - cu1 );
280 T cv0( 1.0 - cv1 );
281 T cw0( 1.0 - cw1 );
282 cu[0] = -0.5*cu1*cu0*cu0; // cubic spline coeffs: u
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; // cubic spline coeffs: v
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; // cubic spline coeffs: w
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 ); // cubic spline grad coeffs: u
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 ); // cubic spline grad coeffs: v
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 ); // cubic spline grad coeffs: w
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;
307 int i, j;
308 for ( j = 0; j < 4; j++ ) {
309 iv = iu;
310 s2 = dv2 = dw2 = 0.0;
311 for ( i = 0; i < 4; i++ ) {
312 iw = iv;
313 s3 = cw[0] * T( map[ iw ] );
314 dw3 = gw[0] * T( map[ iw ] );
315 iw.next_w();
316 s3 += cw[1] * T( map[ iw ] );
317 dw3 += gw[1] * T( map[ iw ] );
318 iw.next_w();
319 s3 += cw[2] * T( map[ iw ] );
320 dw3 += gw[2] * T( map[ iw ] );
321 iw.next_w();
322 s3 += cw[3] * T( map[ iw ] );
323 dw3 += gw[3] * T( map[ iw ] );
324 s2 += cv[i] * s3;
325 dv2 += gv[i] * s3;
326 dw2 += cv[i] * dw3;
327 iv.next_v();
328 }
329 s1 += cu[j] * s2;
330 du1 += gu[j] * s2;
331 dv1 += cu[j] * dv2;
332 dw1 += cu[j] * dw2;
333 iu.next_u();
334 }
335 val = s1;
336 grad = Grad_map<T>( du1, dv1, dw1 );
337 }
338
339
347 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 )
348 {
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 );
360 T cu0( 1.0 - cu1 );
361 T cv0( 1.0 - cv1 );
362 T cw0( 1.0 - cw1 );
363 cu[0] = -0.5*cu1*cu0*cu0; // cubic spline coeffs: u
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; // cubic spline coeffs: v
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; // cubic spline coeffs: w
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 ); // cubic spline grad coeffs: u
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 ); // cubic spline grad coeffs: v
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 ); // cubic spline grad coeffs: w
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; // cubic spline curv coeffs: u
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; // cubic spline curv coeffs: v
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; // cubic spline curv coeffs: w
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;
400 int i, j;
401 for ( j = 0; j < 4; j++ ) {
402 iv = iu;
403 s2 = dv2 = dw2 = dvw2 = dvv2 = dww2 = 0.0;
404 for ( i = 0; i < 4; i++ ) {
405 iw = iv;
406 s3 = cw[0] * T( map[ iw ] );
407 dw3 = gw[0] * T( map[ iw ] );
408 dww3 = ggw[0] * T( map[ iw ] );
409 iw.next_w();
410 s3 += cw[1] * T( map[ iw ] );
411 dw3 += gw[1] * T( map[ iw ] );
412 dww3 += ggw[1] * T( map[ iw ] );
413 iw.next_w();
414 s3 += cw[2] * T( map[ iw ] );
415 dw3 += gw[2] * T( map[ iw ] );
416 dww3 += ggw[2] * T( map[ iw ] );
417 iw.next_w();
418 s3 += cw[3] * T( map[ iw ] );
419 dw3 += gw[3] * T( map[ iw ] );
420 dww3 += ggw[3] * T( map[ iw ] );
421 s2 += cv[i] * s3;
422 dv2 += gv[i] * s3;
423 dw2 += cv[i] * dw3;
424 dvw2 += gv[i] * dw3;
425 dvv2 += ggv[i] * s3;
426 dww2 += cv[i] * dww3;
427 iv.next_v();
428 }
429 s1 += cu[j] * s2;
430 du1 += gu[j] * s2;
431 dv1 += cu[j] * dv2;
432 dw1 += cu[j] * dw2;
433 duv1 += gu[j] * dv2;
434 duw1 += gu[j] * dw2;
435 dvw1 += cu[j] * dvw2;
436 duu1 += ggu[j] * s2;
437 dvv1 += cu[j] * dvv2;
438 dww1 += cu[j] * dww2;
439 iu.next_u();
440 }
441 val = s1;
442 grad = Grad_map<T>( du1, dv1, dw1 );
443 curv = Curv_map<T>( Mat33<T>( duu1, duv1, duw1,
444 duv1, dvv1, dvw1,
445 duw1, dvw1, dww1 ) );
446 }
447
448
449} // namespace clipper
450
451
452#endif
Grid coordinate.
Definition coords.h:237
const int & v() const
get v
Definition coords.h:249
const int & w() const
get w
Definition coords.h:250
const int & u() const
get u
Definition coords.h:248
map coordinate: this is like Coord_grid, but non-integer
Definition coords.h:388
const ftype & u() const
get u
Definition coords.h:408
const ftype & v() const
get v
Definition coords.h:409
Coord_grid coord_grid() const
return integer Coord_grid nearest this coordinate
Definition coords.h:403
const ftype & w() const
get w
Definition coords.h:410
Coord_grid floor() const
return integer Coord_grid below this coordinate
Definition coords.h:405
map coordinate curvatures, with respect to grid u,v,w
Definition derivs.h:146
map coordinate gradient, with respect to grid u,v,w
Definition derivs.h:102
Wrapper class for third-order (cubic) interpolation fns.
Definition map_interp.h:109
static int order()
Order of interpolant.
Definition map_interp.h:115
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
static void interp_grad(const M &map, const Coord_map &pos, T &val, Grad_map< T > &grad)
Definition map_interp.h:267
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
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
Wrapper class for first-order (linear) interpolation fns.
Definition map_interp.h:88
static int order()
Order of interpolant.
Definition map_interp.h:92
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
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
Wrapper class for zeroth-order (nearest neighbour) interpolation fns.
Definition map_interp.h:67
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
static int order()
Order of interpolant.
Definition map_interp.h:71
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
3x3-matrix class
Definition clipper_types.h:183
ftype64 ftype
ftype definition for floating point representation
Definition clipper_precision.h:58