USML
data_grid_bathy.h
1 
6 #ifndef USML_TYPES_DATA_GRID_BATHY_H
7 #define USML_TYPES_DATA_GRID_BATHY_H
8 
9 #include <usml/types/data_grid.h>
10 #include <usml/types/seq_vector.h>
11 
12 //#define FAST_GRID_DEBUG
13 //#define DERV_CONSTRUCT
14 
15 namespace usml {
16 namespace types {
19 
38 class USML_DECLSPEC data_grid_bathy: public data_grid<double, 2> {
39 
40 public:
51  data_grid_bathy(const data_grid<double, 2>& grid, bool copy_data = true) :
52  data_grid<double, 2>(grid, copy_data), _bicubic_coeff(16, 1),
53  _field(16, 1), _xyloc(1, 16), _result_pchip(1, 1), _value(4, 4),
54  _kmin(0u), _k0max(_axis[0]->size() - 1u), _k1max(_axis[1]->size() - 1u)
55  {
56  // Construct the inverse bicubic interpolation coefficient matrix
57  _inv_bicubic_coeff = zero_matrix<double>(16, 16);
58  _inv_bicubic_coeff(0, 0) = 1;
59  _inv_bicubic_coeff(1, 8) = 1;
60  _inv_bicubic_coeff(2, 0) = -3;
61  _inv_bicubic_coeff(2, 1) = 3;
62  _inv_bicubic_coeff(2, 8) = -2;
63  _inv_bicubic_coeff(2, 9) = -1;
64  _inv_bicubic_coeff(3, 0) = 2;
65  _inv_bicubic_coeff(3, 1) = -2;
66  _inv_bicubic_coeff(3, 8) = _inv_bicubic_coeff(3, 9) = 1;
67  _inv_bicubic_coeff(4, 4) = 1;
68  _inv_bicubic_coeff(5, 12) = 1;
69  _inv_bicubic_coeff(6, 4) = -3;
70  _inv_bicubic_coeff(6, 5) = 3;
71  _inv_bicubic_coeff(6, 12) = -2;
72  _inv_bicubic_coeff(6, 13) = -1;
73  _inv_bicubic_coeff(7, 4) = 2;
74  _inv_bicubic_coeff(7, 5) = -2;
75  _inv_bicubic_coeff(7, 12) = _inv_bicubic_coeff(7, 13) = 1;
76  _inv_bicubic_coeff(8, 0) = -3;
77  _inv_bicubic_coeff(8, 2) = 3;
78  _inv_bicubic_coeff(8, 4) = -2;
79  _inv_bicubic_coeff(8, 6) = -1;
80  _inv_bicubic_coeff(9, 8) = -3;
81  _inv_bicubic_coeff(9, 10) = 3;
82  _inv_bicubic_coeff(9, 12) = -2;
83  _inv_bicubic_coeff(9, 14) = -1;
84  _inv_bicubic_coeff(10, 0) = _inv_bicubic_coeff(10, 3) = 9;
85  _inv_bicubic_coeff(10, 1) = _inv_bicubic_coeff(10, 2) = -9;
86  _inv_bicubic_coeff(10, 4) = _inv_bicubic_coeff(10, 8) = 6;
87  _inv_bicubic_coeff(10, 5) = _inv_bicubic_coeff(10, 10) = -6;
88  _inv_bicubic_coeff(10, 6) = _inv_bicubic_coeff(10, 9) = 3;
89  _inv_bicubic_coeff(10, 7) = _inv_bicubic_coeff(10, 11) = -3;
90  _inv_bicubic_coeff(10, 12) = 4;
91  _inv_bicubic_coeff(10, 13) = _inv_bicubic_coeff(10, 14) = 2;
92  _inv_bicubic_coeff(10, 15) = 1;
93  _inv_bicubic_coeff(11, 0) = _inv_bicubic_coeff(11, 3) = -6;
94  _inv_bicubic_coeff(11, 1) = _inv_bicubic_coeff(11, 2) = 6;
95  _inv_bicubic_coeff(11, 6) = _inv_bicubic_coeff(11, 12) =
96  _inv_bicubic_coeff(11, 13) = -2;
97  _inv_bicubic_coeff(11, 4) = -4;
98  _inv_bicubic_coeff(11, 5) = 4;
99  _inv_bicubic_coeff(11, 7) = 2;
100  _inv_bicubic_coeff(11, 8) = _inv_bicubic_coeff(11, 9) = -3;
101  _inv_bicubic_coeff(11, 10) = _inv_bicubic_coeff(11, 11) = 3;
102  _inv_bicubic_coeff(11, 14) = _inv_bicubic_coeff(11, 15) = -1;
103  _inv_bicubic_coeff(12, 0) = 2;
104  _inv_bicubic_coeff(12, 2) = -2;
105  _inv_bicubic_coeff(12, 4) = _inv_bicubic_coeff(12, 6) = 1;
106  _inv_bicubic_coeff(13, 8) = 2;
107  _inv_bicubic_coeff(13, 10) = -2;
108  _inv_bicubic_coeff(13, 12) = _inv_bicubic_coeff(13, 14) = 1;
109  _inv_bicubic_coeff(14, 0) = _inv_bicubic_coeff(14, 3) = -6;
110  _inv_bicubic_coeff(14, 1) = _inv_bicubic_coeff(14, 2) = 6;
111  _inv_bicubic_coeff(14, 4) = _inv_bicubic_coeff(14, 6) = -3;
112  _inv_bicubic_coeff(14, 5) = _inv_bicubic_coeff(14, 7) = 3;
113  _inv_bicubic_coeff(14, 8) = -4;
114  _inv_bicubic_coeff(14, 10) = 4;
115  _inv_bicubic_coeff(14, 9) = _inv_bicubic_coeff(14, 12) =
116  _inv_bicubic_coeff(14, 14) = -2;
117  _inv_bicubic_coeff(14, 11) = 2;
118  _inv_bicubic_coeff(14, 13) = _inv_bicubic_coeff(14, 15) = -1;
119  _inv_bicubic_coeff(15, 0) = _inv_bicubic_coeff(15, 3) = 4;
120  _inv_bicubic_coeff(15, 1) = _inv_bicubic_coeff(15, 2) = -4;
121  _inv_bicubic_coeff(15, 4) = _inv_bicubic_coeff(15, 6) =
122  _inv_bicubic_coeff(15, 8) = _inv_bicubic_coeff(15, 9) = 2;
123  _inv_bicubic_coeff(15, 5) = _inv_bicubic_coeff(15, 7) =
124  _inv_bicubic_coeff(15, 10) = _inv_bicubic_coeff(15, 11) = -2;
125  _inv_bicubic_coeff(15, 12) = _inv_bicubic_coeff(15, 13) =
126  _inv_bicubic_coeff(15, 14) = _inv_bicubic_coeff(15, 15) = 1;
127 
128  _fast_index[0] = 0;
129  _fast_index[1] = 0;
130 
131  //Pre-construct increments for all intervals once to save time
132  matrix<double> inc_x(_k0max + 1u, 1);
133  for (unsigned i = 0; i < _k0max + 1u; ++i) {
134  if (i == 0 || i == _k0max) {
135  inc_x(i, 0) = 2;
136  } else {
137  inc_x(i, 0) = (_axis[0]->increment(i - 1)
138  + _axis[0]->increment(i + 1)) / _axis[0]->increment(i);
139  }
140  }
141  matrix<double> inc_y(_k1max + 1u, 1);
142  for (unsigned i = 0; i < _k1max + 1u; ++i) {
143  if (i == 0 || i == _k1max) {
144  inc_y(i, 0) = 2;
145  } else {
146  inc_y(i, 0) = (_axis[1]->increment(i - 1)
147  + _axis[1]->increment(i + 1)) / _axis[1]->increment(i);
148  }
149  }
150 
151  // Pre-construct all derivatives and cross-dervs once to save time
152  _derv_x = matrix<double>(_k0max + 1u, _k1max + 1u);
153  _derv_y = matrix<double>(_k0max + 1u, _k1max + 1u);
154  _derv_x_y = matrix<double>(_k0max + 1u, _k1max + 1u);
155  for (unsigned i = 0; i < _k0max + 1u; ++i) {
156  for (unsigned j = 0; j < _k1max + 1u; ++j) {
157  if (i < 1 && j < 1) { //top-left corner
158  #ifdef DERV_CONSTRUCT
159  cout << "***Condition: i<1 && j<1***" << endl;
160  #endif
161  _derv_x(i, j) = (data_2d(i + 1, j) - data_2d(i, j))
162  / inc_x(i, 0);
163  _derv_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j))
164  / inc_y(j, 0);
165  _derv_x_y(i, j) = (data_2d(i + 1, j + 1) - data_2d(i + 1, j)
166  - data_2d(i, j + 1) + data_2d(i, j))
167  / (inc_x(i, 0) * inc_y(j, 0));
168  } else if (i == _k0max && j == _k1max) { //bottom-right corner
169  #ifdef DERV_CONSTRUCT
170  cout << "***Condition: i==_k0max && j==_k1max***" << endl;
171  #endif
172  _derv_x(i, j) = (data_2d(i, j) - data_2d(i - 1, j))
173  / inc_x(i, 0);
174  _derv_y(i, j) = (data_2d(i, j) - data_2d(i, j - 1))
175  / inc_y(j, 0);
176  _derv_x_y(i, j) = (data_2d(i, j) - data_2d(i, j - 1)
177  - data_2d(i - 1, j) + data_2d(i - 1, j - 1))
178  / (inc_x(i, 0) * inc_y(j, 0));
179  } else if (i < 1 && j == _k1max) { //top-right corner
180  #ifdef DERV_CONSTRUCT
181  cout << "***Condition: i<1 && j==_k1max***" << endl;
182  #endif
183  _derv_x(i, j) = (data_2d(i + 1, j) - data_2d(i, j))
184  / inc_x(i, 0);
185  _derv_y(i, j) = (data_2d(i, j) - data_2d(i, j - 1))
186  / inc_y(j, 0);
187  _derv_x_y(i, j) = (data_2d(i + 1, j) - data_2d(i + 1, j - 1)
188  - data_2d(i, j) + data_2d(i, j - 1))
189  / (inc_x(i, 0) * inc_y(j, 0));
190  } else if (j < 1 && i == _k0max) { //bottom-left corner
191  #ifdef DERV_CONSTRUCT
192  cout << "***Condition: j<1 && i==_k0max***" << endl;
193  #endif
194  _derv_x(i, j) = (data_2d(i, j) - data_2d(i - 1, j))
195  / inc_x(i, 0);
196  _derv_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j))
197  / inc_y(j, 0);
198  _derv_x_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j)
199  - data_2d(i - 1, j + 1) + data_2d(i - 1, j))
200  / (inc_x(i, 0) * inc_y(j, 0));
201  } else if (i < 1 && (1 <= j && j < _k1max)) { //top row
202  #ifdef DERV_CONSTRUCT
203  cout << "***Condition: i<1 && (1<=j && j<_k1max)***" << endl;
204  #endif
205  _derv_x(i, j) = (data_2d(i + 1, j) - data_2d(i, j))
206  / inc_x(i, 0);
207  _derv_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j - 1))
208  / inc_y(j, 0);
209  _derv_x_y(i, j) = (data_2d(i + 1, j + 1)
210  - data_2d(i + 1, j - 1) - data_2d(i, j + 1)
211  + data_2d(i, j - 1)) / (inc_x(i, 0) * inc_y(j, 0));
212  } else if (j < 1 && (1 <= i && i < _k0max)) { //left most column
213  #ifdef DERV_CONSTRUCT
214  cout << "***Condition: j<1 && (1<=i && i<_k0max)***" << endl;
215  #endif
216  _derv_x(i, j) = (data_2d(i + 1, j) - data_2d(i - 1, j))
217  / inc_x(i, 0);
218  _derv_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j))
219  / inc_y(j, 0);
220  _derv_x_y(i, j) = (data_2d(i + 1, j + 1) - data_2d(i + 1, j)
221  - data_2d(i - 1, j + 1) + data_2d(i - 1, j))
222  / (inc_x(i, 0) * inc_y(j, 0));
223  } else if (j == _k1max && (1 <= i && i < _k0max)) { //right most column
224  #ifdef DERV_CONSTRUCT
225  cout << "***Condition: j>_k1max && (1<=i && i<_k0max)***" << endl;
226  #endif
227  _derv_x(i, j) = (data_2d(i + 1, j) - data_2d(i - 1, j))
228  / inc_x(i, 0);
229  _derv_y(i, j) = (data_2d(i, j) - data_2d(i, j - 1))
230  / inc_y(j, 0);
231  _derv_x_y(i, j) = (data_2d(i + 1, j) - data_2d(i + 1, j - 1)
232  - data_2d(i - 1, j) + data_2d(i - 1, j - 1))
233  / (inc_x(i, 0) * inc_y(j, 0));
234  } else if (i == _k0max && (1 <= j && j < _k1max)) { //bottom row
235  #ifdef DERV_CONSTRUCT
236  cout << "***Condition: i>_k0max && (1<=j && j<_k1max)***" << endl;
237  #endif
238  _derv_x(i, j) = (data_2d(i, j) - data_2d(i - 1, j))
239  / inc_x(i, 0);
240  _derv_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j - 1))
241  / inc_y(j, 0);
242  _derv_x_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j - 1)
243  - data_2d(i - 1, j + 1) + data_2d(i - 1, j - 1))
244  / (inc_x(i, 0) * inc_y(j, 0));
245  } else { //inside, restrictive
246  _derv_x(i, j) = (data_2d(i + 1, j) - data_2d(i - 1, j))
247  / inc_x(i, 0);
248  _derv_y(i, j) = (data_2d(i, j + 1) - data_2d(i, j - 1))
249  / inc_y(j, 0);
250  _derv_x_y(i, j) = (data_2d(i + 1, j + 1)
251  - data_2d(i + 1, j - 1) - data_2d(i - 1, j + 1)
252  + data_2d(i - 1, j - 1))
253  / (inc_x(i, 0) * inc_y(j, 0));
254  }
255  } //end for-loop in j
256  } //end for-loop in i
257  }// end constructor
258 
272  double interpolate(double* location, double* derivative = NULL) {
273 
274  double result = 0;
275  // find the interval index in each dimension
276 
277  for (unsigned dim = 0; dim < 2; ++dim) {
278 
279  // limit interpolation to axis domain if _edge_limit turned on
280 
281  if (edge_limit(dim)) {
282  double a = *(_axis[dim]->begin());
283  double b = *(_axis[dim]->rbegin());
284  double inc = _axis[dim]->increment(0);
285  if (inc < 0) { // a > b
286  if (location[dim] >= a) { //left of the axis
287  location[dim] = a;
288  _offset[dim] = 0;
289  } else if (location[dim] <= b) { //right of the axis
290  location[dim] = b;
291  _offset[dim] = _axis[dim]->size() - 2;
292  } else {
293  _offset[dim] = _axis[dim]->find_index(location[dim]); //somewhere in-between the endpoints of the axis
294  }
295  }
296  if (inc > 0) { // a < b
297  if (location[dim] <= a) { //left of the axis
298  location[dim] = a;
299  _offset[dim] = 0;
300  } else if (location[dim] >= b) { //right of the axis
301  location[dim] = b;
302  _offset[dim] = _axis[dim]->size() - 2;
303  } else {
304  _offset[dim] = _axis[dim]->find_index(location[dim]); //somewhere in-between the endpoints of the axis
305  }
306  }
307 
308  // allow extrapolation if _edge_limit turned off
309 
310  } else {
311  _offset[dim] = _axis[dim]->find_index(location[dim]);
312  }
313  }
314 
315  switch (interp_type(0)) {
316 
317  //****nearest****
318  case -1:
319  for (int dim = 0; dim < 2; ++dim) {
320  double inc = _axis[dim]->increment(0);
321  double u = abs(location[dim] - (*_axis[dim])(_offset[dim]))
322  / inc;
323  if (u < 0.5) {
324  _fast_index[dim] = _offset[dim];
325  } else {
326  _fast_index[dim] = _offset[dim] + 1;
327  }
328  }
329  if (derivative)
330  derivative[0] = derivative[1] = 0;
331  return data(_fast_index);
332  break;
333 
334  //****linear****
335  case 0:
336  double f11, f21, f12, f22, x_diff, y_diff;
337  double x, x1, x2, y, y1, y2;
338 
339  x = location[0];
340  x1 = (*_axis[0])(_offset[0]);
341  x2 = (*_axis[0])(_offset[0] + 1);
342  y = location[1];
343  y1 = (*_axis[1])(_offset[1]);
344  y2 = (*_axis[1])(_offset[1] + 1);
345  f11 = data(_offset);
346  _fast_index[0] = _offset[0] + 1;
347  _fast_index[1] = _offset[1];
348  f21 = data(_fast_index);
349  _fast_index[0] = _offset[0];
350  _fast_index[1] = _offset[1] + 1;
351  f12 = data(_fast_index);
352  _fast_index[0] = _offset[0] + 1;
353  _fast_index[1] = _offset[1] + 1;
354  f22 = data(_fast_index);
355  x_diff = x2 - x1;
356  y_diff = y2 - y1;
357  result = (f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y)
358  + f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1))
359  / (x_diff * y_diff);
360  if (derivative) {
361  derivative[0] = (f21 * (y2 - y) - f11 * (y2 - y)
362  + f22 * (y - y1) - f12 * (y - y1)) / (x_diff * y_diff);
363  derivative[1] = (f12 * (x2 - x) - f11 * (x2 - x)
364  + f22 * (x - x1) - f21 * (x - x1)) / (x_diff * y_diff);
365  }
366  return result;
367  break;
368 
369  //****pchip****
370  case 1:
371  result = fast_pchip(_offset, location, derivative);
372  return result;
373  break;
374 
375  default:
376  throw std::invalid_argument(
377  "Interp must be NEAREST, LINEAR, or PCHIP");
378  break;
379  }
380  } // end interpolate at a single location
381 
395  void interpolate(const matrix<double>& x, const matrix<double>& y,
396  matrix<double>* result, matrix<double>* dx = NULL,
397  matrix<double>* dy = NULL) {
398  double location[2];
399  double derivative[2];
400  for (unsigned n = 0; n < x.size1(); ++n) {
401  for (unsigned m = 0; m < x.size2(); ++m) {
402  location[0] = x(n, m);
403  location[1] = y(n, m);
404  if (dx == NULL || dy == NULL) {
405  (*result)(n, m) = (double) interpolate(location);
406  } else {
407  (*result)(n, m) = (double) interpolate(location,
408  derivative);
409  (*dx)(n, m) = (double) derivative[0];
410  (*dy)(n, m) = (double) derivative[1];
411  }
412  }
413  }
414  } // end Interpolate at a series of locations.
415 
416 private:
417 
419  inline double data_2d(unsigned row, unsigned col) {
420  unsigned grid_index[2];
421  grid_index[0] = row;
422  grid_index[1] = col;
423  return data(grid_index);
424  }
425 
478  double fast_pchip(const unsigned* interp_index, double* location,
479  double* derivative = NULL) {
480  int k0 = interp_index[0];
481  int k1 = interp_index[1];
482  double norm0, norm1;
483 
484  // Checks for boundaries of the axes
485  norm0 = (*_axis[0])(k0 + 1) - (*_axis[0])(k0);
486  norm1 = (*_axis[1])(k1 + 1) - (*_axis[1])(k1);
487  for (int i = -1; i < 3; ++i) {
488  for (int j = -1; j < 3; ++j) {
489  //get appropriate data when at boundaries
490  if ((k0 + i) >= _k0max) {
491  _fast_index[0] = _k0max;
492  } else if ((k0 + i) <= _kmin) {
493  _fast_index[0] = _kmin;
494  } else {
495  _fast_index[0] = k0 + i;
496  }
497  //get appropriate data when at boundaries
498  if ((k1 + j) >= _k1max) {
499  _fast_index[1] = _k1max;
500  } else if ((k1 + j) <= _kmin) {
501  _fast_index[1] = _kmin;
502  } else {
503  _fast_index[1] = k1 + j;
504  }
505  _value(i + 1, j + 1) = data(_fast_index);
506  } //end for-loop in j
507  } //end for-loop in i
508 
509 #ifdef FAST_GRID_DEBUG
510  cout << "loc0: " << location[0] << " loc1: " << location[1] << endl;
511  cout << "offset0: " << k0 << " offset1: " << k1 << endl;
512  cout << "axis0: " << (*_axis[0])(k0) << " axis1: " << (*_axis[1])(k1) << endl;
513  cout << "data value at offset: " << ( (data(interp_index) > 1e6) ?
514  data(interp_index) - wposition::earth_radius : data(interp_index) ) << endl;
515  cout << "_value: " << _value << endl;
516  cout << "_derv_x: " << _derv_x << endl;
517  cout << "_derv_y: " << _derv_y << endl;
518  cout << "_derv_x_y: " << _derv_x_y << endl;
519 #endif
520 
521  // Construct the _field matrix
522  _field(0, 0) = _value(1, 1); //f(0,0)
523  _field(1, 0) = _value(1, 2); //f(0,1)
524  _field(2, 0) = _value(2, 1); //f(1,0)
525  _field(3, 0) = _value(2, 2); //f(1,1)
526  _field(4, 0) = _derv_x(k0, k1); //f_x(0,0)
527  _field(5, 0) = _derv_x(k0, k1 + 1); //f_x(0,1)
528  _field(6, 0) = _derv_x(k0 + 1, k1); //f_x(1,0)
529  _field(7, 0) = _derv_x(k0 + 1, k1 + 1); //f_x(1,1)
530  _field(8, 0) = _derv_y(k0, k1); //f_y(0,0)
531  _field(9, 0) = _derv_y(k0, k1 + 1); //f_y(0,1)
532  _field(10, 0) = _derv_y(k0 + 1, k1); //f_y(1,0)
533  _field(11, 0) = _derv_y(k0 + 1, k1 + 1); //f_y(1,1)
534  _field(12, 0) = _derv_x_y(k0, k1); //f_x_y(0,0)
535  _field(13, 0) = _derv_x_y(k0, k1 + 1); //f_x_y(0,1)
536  _field(14, 0) = _derv_x_y(k0 + 1, k1); //f_x_y(1,0)
537  _field(15, 0) = _derv_x_y(k0 + 1, k1 + 1); //f_x_y(1,1)
538 
539  // Construct the coefficients of the bicubic interpolation
540  _bicubic_coeff = prod(_inv_bicubic_coeff, _field);
541 
542  // Create the power series of the interpolation formula before hand for speed
543  double x_inv = location[0] - (*_axis[0])(k0);
544  double y_inv = location[1] - (*_axis[1])(k1);
545 
546 #ifdef FAST_GRID_DEBUG
547  cout << "_field: " << _field << endl;
548  cout << "_bicubic_coeff: " << _bicubic_coeff << endl;
549  cout << "x_inv/norm0: " << x_inv/norm0 << "\ty_inv/norm1: " << y_inv/norm1 << endl;
550 #endif
551 
552  _xyloc(0, 0) = 1;
553  _xyloc(0, 1) = y_inv / norm1;
554  _xyloc(0, 2) = _xyloc(0, 1) * _xyloc(0, 1);
555  _xyloc(0, 3) = _xyloc(0, 2) * _xyloc(0, 1);
556  _xyloc(0, 4) = x_inv / norm0;
557  _xyloc(0, 5) = _xyloc(0, 4) * _xyloc(0, 1);
558  _xyloc(0, 6) = _xyloc(0, 4) * _xyloc(0, 2);
559  _xyloc(0, 7) = _xyloc(0, 4) * _xyloc(0, 3);
560  _xyloc(0, 8) = _xyloc(0, 4) * _xyloc(0, 4);
561  _xyloc(0, 9) = _xyloc(0, 8) * _xyloc(0, 1);
562  _xyloc(0, 10) = _xyloc(0, 8) * _xyloc(0, 2);
563  _xyloc(0, 11) = _xyloc(0, 8) * _xyloc(0, 3);
564  _xyloc(0, 12) = _xyloc(0, 8) * _xyloc(0, 4);
565  _xyloc(0, 13) = _xyloc(0, 12) * _xyloc(0, 1);
566  _xyloc(0, 14) = _xyloc(0, 12) * _xyloc(0, 2);
567  _xyloc(0, 15) = _xyloc(0, 12) * _xyloc(0, 3);
568 
569  _result_pchip = prod(_xyloc, _bicubic_coeff);
570  if (derivative) {
571  derivative[0] = 0;
572  derivative[1] = 0;
573  for (int i = 1; i < 4; ++i) {
574  for (int j = 0; j < 4; ++j) {
575  derivative[0] += i * _bicubic_coeff(i * 4 + j, 0)
576  * _xyloc(0, 4 * (i - 1)) * _xyloc(0, j);
577  }
578  }
579  for (int i = 0; i < 4; ++i) {
580  for (int j = 1; j < 4; ++j) {
581  derivative[1] += j * _bicubic_coeff(i * 4 + j, 0)
582  * _xyloc(0, 4 * i) * _xyloc(0, j - 1);
583  }
584  }
585  }
586  return _result_pchip(0, 0);
587  }
588 
589  //***********************************************************/
590  // Private data members
591 
597  c_matrix<double, 16, 16> _inv_bicubic_coeff;
598 
606  c_matrix<double, 16, 1> _bicubic_coeff;
607  c_matrix<double, 16, 1> _field;
608  c_matrix<double, 1, 16> _xyloc;
609  c_matrix<double, 1, 1> _result_pchip;
610  c_matrix<double, 4, 4> _value;
611  matrix<double> _derv_x;
612  matrix<double> _derv_y;
613  matrix<double> _derv_x_y;
614  unsigned _fast_index[2];
615  const int _kmin;
616  const int _k0max;
617  const int _k1max;
618 
619 }; // end data_grid_bathy
620 
621 } // end of namespace types
622 } // end of namespace usml
623 
624 #endif
c_matrix< double, 1, 1 > _result_pchip
Definition: data_grid_bathy.h:609
static double earth_radius
Local radius of curvature in the area of operations.
Definition: wposition.h:100
const int _kmin
Definition: data_grid_bathy.h:615
matrix< double > _derv_y
Definition: data_grid_bathy.h:612
c_matrix< double, 16, 1 > _bicubic_coeff
Create variables that are used multiple times through out the course of multiple calls per instance...
Definition: data_grid_bathy.h:606
BOOST_UBLAS_INLINE matrix_unary1_traits< E, scalar_abs< typename E::value_type > >::result_type abs(const matrix_expression< E > &e)
Magnitude of a complex matrix.
Definition: matrix_math.h:257
Implements fast calculations for data_grids using a non-recursive engine on interpolation.
Definition: data_grid_bathy.h:38
const int _k0max
Definition: data_grid_bathy.h:616
const int _k1max
Definition: data_grid_bathy.h:617
void interpolate(const matrix< double > &x, const matrix< double > &y, matrix< double > *result, matrix< double > *dx=NULL, matrix< double > *dy=NULL)
Overrides the interpolate function within data_grid using the non-recursive formula.
Definition: data_grid_bathy.h:395
matrix< double > _derv_x_y
Definition: data_grid_bathy.h:613
data_grid_bathy(const data_grid< double, 2 > &grid, bool copy_data=true)
Constructor - Creates a fast interpolation grid from an existing data_grid.
Definition: data_grid_bathy.h:51
c_matrix< double, 16, 16 > _inv_bicubic_coeff
Matrix that is the inverse of the bicubic coefficients This matrix will be used to construct the bicu...
Definition: data_grid_bathy.h:597
N-dimensional data set and its associated axes.
Definition: data_grid.h:64
double interpolate(double *location, double *derivative=NULL)
Overrides the interpolate function within data_grid using the non-recursive formula.
Definition: data_grid_bathy.h:272
double data_2d(unsigned row, unsigned col)
Utility accessor function for data grid values.
Definition: data_grid_bathy.h:419
c_matrix< double, 1, 16 > _xyloc
Definition: data_grid_bathy.h:608
double fast_pchip(const unsigned *interp_index, double *location, double *derivative=NULL)
A non-recursive version of the Piecewise Cubic Hermite polynomial (PCHIP) specific to the 2-dimension...
Definition: data_grid_bathy.h:478
matrix< double > _derv_x
Definition: data_grid_bathy.h:611
c_matrix< double, 4, 4 > _value
Definition: data_grid_bathy.h:610
c_matrix< double, 16, 1 > _field
Definition: data_grid_bathy.h:607