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 
51 namespace 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
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