GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
InterpSpline.c
Go to the documentation of this file.
1 
2 /***********************************************************************
3  *
4  * MODULE: lidarlib
5  *
6  * AUTHOR(S): Roberto Antolin
7  *
8  * PURPOSE: LIDAR Spline Interpolation
9  *
10  * COPYRIGHT: (C) 2006 by Politecnico di Milano -
11  * Polo Regionale di Como
12  *
13  * This program is free software under the
14  * GNU General Public License (>=v2).
15  * Read the file COPYING that comes with GRASS
16  * for details.
17  *
18  **************************************************************************/
19 
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <float.h>
23 #include <math.h>
24 #include <string.h>
25 #include <grass/lidar.h>
26 
27 /*----------------------------------------------------------------------------*/
28 /* Abscissa node index computation */
29 
30 void node_x(double x, int *i_x, double *csi_x, double xMin, double deltaX)
31 {
32 
33  *i_x = (int)((x - xMin) / deltaX);
34  *csi_x = (double)((x - xMin) - (*i_x * deltaX));
35 
36  return;
37 }
38 
39 /*----------------------------------------------------------------------------*/
40 /* Ordinate node index computation */
41 
42 void node_y(double y, int *i_y, double *csi_y, double yMin, double deltaY)
43 {
44 
45  *i_y = (int)((y - yMin) / deltaY);
46  *csi_y = (double)((y - yMin) - (*i_y * deltaY));
47 
48  return;
49 }
50 
51 /*----------------------------------------------------------------------------*/
52 /* Node order computation */
53 
54 int order(int i_x, int i_y, int yNum)
55 {
56 
57  return (i_y + i_x * yNum);
58 }
59 
60 /*----------------------------------------------------------------------------*/
61 /* Design matrix coefficients computation */
62 
63 double phi_3(double csi)
64 {
65 
66  return ((pow(2 - csi, 3.) - pow(1 - csi, 3.) * 4) / 6.);
67 }
68 
69 double phi_4(double csi)
70 {
71 
72  return (pow(2 - csi, 3.) / 6.);
73 }
74 
75 double phi_33(double csi_x, double csi_y)
76 {
77 
78  return (phi_3(csi_x) * phi_3(csi_y));
79 }
80 
81 double phi_34(double csi_x, double csi_y)
82 {
83 
84  return (phi_3(csi_x) * phi_4(csi_y));
85 }
86 
87 double phi_43(double csi_x, double csi_y)
88 {
89 
90  return (phi_4(csi_x) * phi_3(csi_y));
91 }
92 
93 double phi_44(double csi_x, double csi_y)
94 {
95 
96  return (phi_4(csi_x) * phi_4(csi_y));
97 }
98 
99 double phi(double csi_x, double csi_y)
100 {
101 
102  return ((1 - csi_x) * (1 - csi_y));
103 }
104 
105 /*----------------------------------------------------------------------------*/
106 /* Normal system computation for bicubic spline interpolation */
107 
108 void normalDefBicubic(double **N, double *TN, double *Q, double **obsVect,
109  double deltaX, double deltaY, int xNum, int yNum,
110  double xMin, double yMin, int obsNum, int parNum,
111  int BW)
112 {
113 
114  int i, k, h, m, n, n0; /* counters */
115  double alpha[4][4]; /* coefficients */
116 
117  int i_x; /* x = (xMin + (i_x * deltaX) + csi_x) */
118  double csi_x;
119 
120  int i_y; /* y = (yMin + (i_y * deltaY) + csi_y) */
121  double csi_y;
122 
123  /*--------------------------------------*/
124  for (k = 0; k < parNum; k++) {
125  for (h = 0; h < BW; h++)
126  N[k][h] = 0.; /* Normal matrix inizialization */
127  TN[k] = 0.; /* Normal vector inizialization */
128  }
129 
130  /*--------------------------------------*/
131 
132  for (i = 0; i < obsNum; i++) {
133 
134  node_x(obsVect[i][0], &i_x, &csi_x, xMin, deltaX);
135  node_y(obsVect[i][1], &i_y, &csi_y, yMin, deltaY);
136 
137  if ((i_x >= -2) && (i_x <= xNum) && (i_y >= -2) && (i_y <= yNum)) {
138 
139  csi_x = csi_x / deltaX;
140  csi_y = csi_y / deltaY;
141 
142  alpha[0][0] = phi_44(1 + csi_x, 1 + csi_y);
143  alpha[0][1] = phi_43(1 + csi_x, csi_y);
144  alpha[0][2] = phi_43(1 + csi_x, 1 - csi_y);
145  alpha[0][3] = phi_44(1 + csi_x, 2 - csi_y);
146 
147  alpha[1][0] = phi_34(csi_x, 1 + csi_y);
148  alpha[1][1] = phi_33(csi_x, csi_y);
149  alpha[1][2] = phi_33(csi_x, 1 - csi_y);
150  alpha[1][3] = phi_34(csi_x, 2 - csi_y);
151 
152  alpha[2][0] = phi_34(1 - csi_x, 1 + csi_y);
153  alpha[2][1] = phi_33(1 - csi_x, csi_y);
154  alpha[2][2] = phi_33(1 - csi_x, 1 - csi_y);
155  alpha[2][3] = phi_34(1 - csi_x, 2 - csi_y);
156 
157  alpha[3][0] = phi_44(2 - csi_x, 1 + csi_y);
158  alpha[3][1] = phi_43(2 - csi_x, csi_y);
159  alpha[3][2] = phi_43(2 - csi_x, 1 - csi_y);
160  alpha[3][3] = phi_44(2 - csi_x, 2 - csi_y);
161 
162  for (k = -1; k <= 2; k++) {
163  for (h = -1; h <= 2; h++) {
164 
165  if (((i_x + k) >= 0) && ((i_x + k) < xNum) &&
166  ((i_y + h) >= 0) && ((i_y + h) < yNum)) {
167  for (m = k; m <= 2; m++) {
168 
169  if (m == k)
170  n0 = h;
171  else
172  n0 = -1;
173 
174  for (n = n0; n <= 2; n++) {
175  if (((i_x + m) >= 0) && ((i_x + m) < xNum) &&
176  ((i_y + n) >= 0) && ((i_y + n) < yNum)) {
177  N[order(i_x + k, i_y + h, yNum)][order
178  (i_x + m,
179  i_y + n,
180  yNum) -
181  order(i_x
182  +
183  k,
184  i_y
185  +
186  h,
187  yNum)]
188  +=
189  alpha[k + 1][h +
190  1] * (1 / Q[i]) *
191  alpha[m + 1][n + 1];
192  /* 1/Q[i] only refers to the variances */
193  }
194  }
195  }
196  TN[order(i_x + k, i_y + h, yNum)] +=
197  obsVect[i][2] * (1 / Q[i]) * alpha[k + 1][h + 1];
198  }
199  }
200  }
201  }
202  }
203 
204  return;
205 }
206 
207 /*----------------------------------------------------------------------------*/
208 /* Normal system correction - Introduzione della correzione dovuta alle
209  pseudosservazioni (Tykonov) - LAPALCIANO - */
210 
211 void nCorrectLapl(double **N, double lambda, int xNum, int yNum,
212  double deltaX, double deltaY)
213 {
214 
215  int i_x, i_y; /* counters */
216  int k, h, m, n, n0; /* counters */
217 
218  double alpha[5][5]; /* coefficients */
219 
220  double lambdaX, lambdaY;
221 
222  /*--------------------------------------*/
223  lambdaX = lambda * (deltaY / deltaX);
224  lambdaY = lambda * (deltaX / deltaY);
225 
226  alpha[0][0] = 0;
227  alpha[0][1] = lambdaX * (1 / 36.); /* There is lambda because Q^(-1) contains 1/(1/lambda) */
228  alpha[0][2] = lambdaX * (1 / 9.);
229  alpha[0][3] = lambdaX * (1 / 36.);
230  alpha[0][4] = 0;
231 
232  alpha[1][0] = lambdaY * (1 / 36.);
233  alpha[1][1] = lambdaX * (1 / 18.) + lambdaY * (1 / 18.);
234  alpha[1][2] = lambdaX * (2 / 9.) - lambdaY * (1 / 6.);
235  alpha[1][3] = lambdaX * (1 / 18.) + lambdaY * (1 / 18.);
236  alpha[1][4] = lambdaY * (1 / 36.);
237 
238  alpha[2][0] = lambdaY * (1 / 9.);
239  alpha[2][1] = -lambdaX * (1 / 6.) + lambdaY * (2 / 9.);
240  alpha[2][2] = -lambdaX * (2 / 3.) - lambdaY * (2 / 3.);
241  alpha[2][3] = -lambdaX * (1 / 6.) + lambdaY * (2 / 9.);
242  alpha[2][4] = lambdaY * (1 / 9.);
243 
244  alpha[3][0] = lambdaY * (1 / 36.);
245  alpha[3][1] = lambdaX * (1 / 18.) + lambdaY * (1 / 18.);
246  alpha[3][2] = lambdaX * (2 / 9.) - lambdaY * (1 / 6.);
247  alpha[3][3] = lambdaX * (1 / 18.) + lambdaY * (1 / 18.);
248  alpha[3][4] = lambdaY * (1 / 36.);
249 
250  alpha[4][0] = 0;
251  alpha[4][1] = lambdaX * (1 / 36.);
252  alpha[4][2] = lambdaX * (1 / 9.);
253  alpha[4][3] = lambdaX * (1 / 36.);
254  alpha[4][4] = 0;
255 
256  for (i_x = 0; i_x < xNum; i_x++) {
257  for (i_y = 0; i_y < yNum; i_y++) {
258 
259  for (k = -2; k <= 2; k++) {
260  for (h = -2; h <= 2; h++) {
261 
262  if (((i_x + k) >= 0) && ((i_x + k) < xNum) &&
263  ((i_y + h) >= 0) && ((i_y + h) < yNum)) {
264 
265  for (m = k; m <= 2; m++) {
266 
267  if (m == k)
268  n0 = h;
269  else
270  n0 = -2;
271 
272  for (n = n0; n <= 2; n++) {
273  if (((i_x + m) >= 0) &&
274  ((i_x + m) <= (xNum - 1)) &&
275  ((i_y + n) >= 0) &&
276  ((i_y + n) <= (yNum - 1))) {
277 
278  if ((alpha[k + 2][h + 2] != 0) &&
279  (alpha[m + 2][n + 2] != 0)) {
280  N[order(i_x + k, i_y + h, yNum)][order
281  (i_x
282  + m,
283  i_y
284  + n,
285  yNum)
286  -
287  order
288  (i_x
289  + k,
290  i_y
291  + h,
292  yNum)]
293  +=
294  alpha[k + 2][h + 2] * alpha[m +
295  2][n +
296  2];
297  }
298  }
299  }
300  }
301  }
302  }
303  }
304  }
305  }
306 
307  return;
308 }
309 
310 /*----------------------------------------------------------------------------*/
311 /* Normal system computation for bilinear spline interpolation */
312 
313 void normalDefBilin(double **N, double *TN, double *Q, double **obsVect,
314  double deltaX, double deltaY, int xNum, int yNum,
315  double xMin, double yMin, int obsNum, int parNum, int BW)
316 {
317 
318  int i, k, h, m, n, n0; /* counters */
319  double alpha[2][2]; /* coefficients */
320 
321  int i_x; /* x = (xMin + (i_x * deltaX) + csi_x) */
322  double csi_x;
323 
324  int i_y; /* y = (yMin + (i_y * deltaY) + csi_y) */
325  double csi_y;
326 
327  /*--------------------------------------*/
328  for (k = 0; k < parNum; k++) {
329  for (h = 0; h < BW; h++)
330  N[k][h] = 0.; /* Normal matrix inizialization */
331  TN[k] = 0.; /* Normal vector inizialization */
332  }
333 
334  /*--------------------------------------*/
335 
336  for (i = 0; i < obsNum; i++) {
337 
338  node_x(obsVect[i][0], &i_x, &csi_x, xMin, deltaX);
339  node_y(obsVect[i][1], &i_y, &csi_y, yMin, deltaY);
340 
341  if ((i_x >= -1) && (i_x < xNum) && (i_y >= -1) && (i_y < yNum)) {
342 
343  csi_x = csi_x / deltaX;
344  csi_y = csi_y / deltaY;
345 
346  alpha[0][0] = phi(csi_x, csi_y);
347  alpha[0][1] = phi(csi_x, 1 - csi_y);
348  alpha[1][0] = phi(1 - csi_x, csi_y);
349  alpha[1][1] = phi(1 - csi_x, 1 - csi_y);
350 
351  for (k = 0; k <= 1; k++) {
352  for (h = 0; h <= 1; h++) {
353 
354  if (((i_x + k) >= 0) && ((i_x + k) <= (xNum - 1)) &&
355  ((i_y + h) >= 0) && ((i_y + h) <= (yNum - 1))) {
356 
357  for (m = k; m <= 1; m++) {
358  if (m == k)
359  n0 = h;
360  else
361  n0 = 0;
362 
363  for (n = n0; n <= 1; n++) {
364  if (((i_x + m) >= 0) && ((i_x + m) < xNum) &&
365  ((i_y + n) >= 0) && ((i_y + n) < yNum)) {
366  N[order(i_x + k, i_y + h, yNum)][order
367  (i_x + m,
368  i_y + n,
369  yNum) -
370  order(i_x
371  +
372  k,
373  i_y
374  +
375  h,
376  yNum)]
377  +=
378  alpha[k][h] * (1 / Q[i]) *
379  alpha[m][n];
380  /* 1/Q[i] only refers to the variances */
381  }
382  }
383  }
384  TN[order(i_x + k, i_y + h, yNum)] +=
385  obsVect[i][2] * (1 / Q[i]) * alpha[k][h];
386  }
387  }
388  }
389  }
390  }
391 
392  return;
393 }
394 
395 
396 /*----------------------------------------------------------------------------*/
397 /* Normal system correction - Introduzione della correzione dovuta alle
398  pseudosservazioni (Tykonov) - GRADIENTE - */
399 #ifdef notdef
400 void nCorrectGrad(double **N, double lambda, int xNum, int yNum,
401  double deltaX, double deltaY)
402 {
403 
404  int i;
405  int parNum;
406 
407  double alpha[3];
408  double lambdaX, lambdaY;
409 
410  lambdaX = lambda * (deltaY / deltaX);
411  lambdaY = lambda * (deltaX / deltaY);
412 
413  parNum = xNum * yNum;
414 
415  alpha[0] = lambdaX / 2. + lambdaY / 2.;
416  alpha[1] = -lambdaX / 4.;
417  alpha[2] = -lambdaY / 4.;
418 
419  for (i = 0; i < parNum; i++) {
420  N[i][0] += alpha[0];
421 
422  if ((i + 2) < parNum)
423  N[i][2] += alpha[2];
424 
425  if ((i + 2 * yNum) < parNum)
426  N[i][2 * yNum] += alpha[1];
427  }
428 }
429 #endif
430 
431 /*1-DELTA discretization */
432 void nCorrectGrad(double **N, double lambda, int xNum, int yNum,
433  double deltaX, double deltaY)
434 {
435 
436  int i;
437  int parNum;
438 
439  double alpha[3];
440  double lambdaX, lambdaY;
441 
442  lambdaX = lambda * (deltaY / deltaX);
443  lambdaY = lambda * (deltaX / deltaY);
444 
445  parNum = xNum * yNum;
446 
447  alpha[0] = 2 * lambdaX + 2 * lambdaY;
448  alpha[1] = -lambdaX;
449  alpha[2] = -lambdaY;
450 
451  for (i = 0; i < parNum; i++) {
452  N[i][0] += alpha[0];
453 
454  if ((i + 1) < parNum)
455  N[i][1] += alpha[2];
456 
457  if ((i + 1 * yNum) < parNum)
458  N[i][1 * yNum] += alpha[1];
459  }
460 
461  return;
462 }
463 
464 /*----------------------------------------------------------------------------*/
465 /* Observations estimation */
466 
467 void obsEstimateBicubic(double **obsV, double *obsE, double *parV,
468  double deltX, double deltY, int xNm, int yNm,
469  double xMi, double yMi, int obsN)
470 {
471 
472  int i, k, h; /* counters */
473  double alpha[4][4]; /* coefficients */
474 
475  int i_x; /* x = (xMin + (i_x * deltaX) + csi_x) */
476  double csi_x;
477 
478  int i_y; /* y = (yMin + (i_y * deltaY) + csi_y) */
479  double csi_y;
480 
481  for (i = 0; i < obsN; i++) {
482 
483  obsE[i] = 0;
484 
485  node_x(obsV[i][0], &i_x, &csi_x, xMi, deltX);
486  node_y(obsV[i][1], &i_y, &csi_y, yMi, deltY);
487 
488  if ((i_x >= -2) && (i_x <= xNm) && (i_y >= -2) && (i_y <= yNm)) {
489 
490  csi_x = csi_x / deltX;
491  csi_y = csi_y / deltY;
492 
493  alpha[0][0] = phi_44(1 + csi_x, 1 + csi_y);
494  alpha[0][1] = phi_43(1 + csi_x, csi_y);
495  alpha[0][2] = phi_43(1 + csi_x, 1 - csi_y);
496  alpha[0][3] = phi_44(1 + csi_x, 2 - csi_y);
497 
498  alpha[1][0] = phi_34(csi_x, 1 + csi_y);
499  alpha[1][1] = phi_33(csi_x, csi_y);
500  alpha[1][2] = phi_33(csi_x, 1 - csi_y);
501  alpha[1][3] = phi_34(csi_x, 2 - csi_y);
502 
503  alpha[2][0] = phi_34(1 - csi_x, 1 + csi_y);
504  alpha[2][1] = phi_33(1 - csi_x, csi_y);
505  alpha[2][2] = phi_33(1 - csi_x, 1 - csi_y);
506  alpha[2][3] = phi_34(1 - csi_x, 2 - csi_y);
507 
508  alpha[3][0] = phi_44(2 - csi_x, 1 + csi_y);
509  alpha[3][1] = phi_43(2 - csi_x, csi_y);
510  alpha[3][2] = phi_43(2 - csi_x, 1 - csi_y);
511  alpha[3][3] = phi_44(2 - csi_x, 2 - csi_y);
512 
513  for (k = -1; k <= 2; k++) {
514  for (h = -1; h <= 2; h++) {
515  if (((i_x + k) >= 0) && ((i_x + k) < xNm) &&
516  ((i_y + h) >= 0) && ((i_y + h) < yNm))
517  obsE[i] +=
518  parV[order(i_x + k, i_y + h, yNm)] * alpha[k +
519  1][h +
520  1];
521  }
522  }
523  }
524  }
525 
526  return;
527 }
528 
529 
530 /*--------------------------------------------------------------------------------------*/
531 /* Data interpolation in a generic point */
532 
533 double dataInterpolateBicubic(double x, double y, double deltaX,
534  double deltaY, int xNum, int yNum, double xMin,
535  double yMin, double *parVect)
536 {
537 
538  double z; /* abscissa, ordinate and associated value */
539 
540  int k, h; /* counters */
541  double alpha[4][4]; /* coefficients */
542 
543  int i_x, i_y; /* x = (xMin + (i_x * deltaX) + csi_x) */
544  double csi_x, csi_y; /* y = (yMin + (i_y * deltaY) + csi_y) */
545 
546 
547  z = 0;
548 
549  node_x(x, &i_x, &csi_x, xMin, deltaX);
550  node_y(y, &i_y, &csi_y, yMin, deltaY);
551 
552  if ((i_x >= -2) && (i_x <= xNum) && (i_y >= -2) && (i_y <= yNum)) {
553 
554  csi_x = csi_x / deltaX;
555  csi_y = csi_y / deltaY;
556 
557  alpha[0][0] = phi_44(1 + csi_x, 1 + csi_y);
558  alpha[0][1] = phi_43(1 + csi_x, csi_y);
559  alpha[0][2] = phi_43(1 + csi_x, 1 - csi_y);
560  alpha[0][3] = phi_44(1 + csi_x, 2 - csi_y);
561 
562  alpha[1][0] = phi_34(csi_x, 1 + csi_y);
563  alpha[1][1] = phi_33(csi_x, csi_y);
564  alpha[1][2] = phi_33(csi_x, 1 - csi_y);
565  alpha[1][3] = phi_34(csi_x, 2 - csi_y);
566 
567  alpha[2][0] = phi_34(1 - csi_x, 1 + csi_y);
568  alpha[2][1] = phi_33(1 - csi_x, csi_y);
569  alpha[2][2] = phi_33(1 - csi_x, 1 - csi_y);
570  alpha[2][3] = phi_34(1 - csi_x, 2 - csi_y);
571 
572  alpha[3][0] = phi_44(2 - csi_x, 1 + csi_y);
573  alpha[3][1] = phi_43(2 - csi_x, csi_y);
574  alpha[3][2] = phi_43(2 - csi_x, 1 - csi_y);
575  alpha[3][3] = phi_44(2 - csi_x, 2 - csi_y);
576 
577  for (k = -1; k <= 2; k++) {
578  for (h = -1; h <= 2; h++) {
579  if (((i_x + k) >= 0) && ((i_x + k) < xNum) && ((i_y + h) >= 0)
580  && ((i_y + h) < yNum))
581  z += parVect[order(i_x + k, i_y + h, yNum)] * alpha[k +
582  1][h +
583  1];
584  }
585  }
586  }
587 
588  return z;
589 }
590 
591 
592 /*----------------------------------------------------------------------------*/
593 /* Observations estimation */
594 
595 void obsEstimateBilin(double **obsV, double *obsE, double *parV, double deltX,
596  double deltY, int xNm, int yNm, double xMi, double yMi,
597  int obsN)
598 {
599 
600  int i, k, h; /* counters */
601  double alpha[2][2]; /* coefficients */
602 
603  int i_x; /* x = (xMin + (i_x * deltaX) + csi_x) */
604  double csi_x;
605 
606  int i_y; /* y = (yMin + (i_y * deltaY) + csi_y) */
607  double csi_y;
608 
609  for (i = 0; i < obsN; i++) {
610 
611  obsE[i] = 0;
612 
613  node_x(obsV[i][0], &i_x, &csi_x, xMi, deltX);
614  node_y(obsV[i][1], &i_y, &csi_y, yMi, deltY);
615 
616  if ((i_x >= -1) && (i_x < xNm) && (i_y >= -1) && (i_y < yNm)) {
617 
618  csi_x = csi_x / deltX;
619  csi_y = csi_y / deltY;
620 
621  alpha[0][0] = phi(csi_x, csi_y);
622  alpha[0][1] = phi(csi_x, 1 - csi_y);
623  alpha[1][0] = phi(1 - csi_x, csi_y);
624  alpha[1][1] = phi(1 - csi_x, 1 - csi_y);
625 
626  for (k = 0; k <= 1; k++) {
627  for (h = 0; h <= 1; h++) {
628  if (((i_x + k) >= 0) && ((i_x + k) < xNm) &&
629  ((i_y + h) >= 0) && ((i_y + h) < yNm))
630  obsE[i] +=
631  parV[order(i_x + k, i_y + h, yNm)] * alpha[k][h];
632  }
633  }
634  }
635  }
636 
637  return;
638 }
639 
640 
641 /*--------------------------------------------------------------------------------------*/
642 /* Data interpolation in a generic point */
643 
644 double dataInterpolateBilin(double x, double y, double deltaX, double deltaY,
645  int xNum, int yNum, double xMin, double yMin,
646  double *parVect)
647 {
648 
649  double z; /* abscissa, ordinate and associated value */
650 
651  int k, h; /* counters */
652  double alpha[2][2]; /* coefficients */
653 
654  int i_x, i_y; /* x = (xMin + (i_x * deltaX) + csi_x) */
655  double csi_x, csi_y; /* y = (yMin + (i_y * deltaY) + csi_y) */
656 
657  z = 0;
658 
659  node_x(x, &i_x, &csi_x, xMin, deltaX);
660  node_y(y, &i_y, &csi_y, yMin, deltaY);
661 
662  if ((i_x >= -1) && (i_x < xNum) && (i_y >= -1) && (i_y < yNum)) {
663 
664  csi_x = csi_x / deltaX;
665  csi_y = csi_y / deltaY;
666 
667  alpha[0][0] = phi(csi_x, csi_y);
668  alpha[0][1] = phi(csi_x, 1 - csi_y);
669 
670  alpha[1][0] = phi(1 - csi_x, csi_y);
671  alpha[1][1] = phi(1 - csi_x, 1 - csi_y);
672 
673  for (k = 0; k <= 1; k++) {
674  for (h = 0; h <= 1; h++) {
675  if (((i_x + k) >= 0) && ((i_x + k) < xNum) && ((i_y + h) >= 0)
676  && ((i_y + h) < yNum))
677  z += parVect[order(i_x + k, i_y + h, yNum)] * alpha[k][h];
678  }
679  }
680  }
681 
682  return z;
683 }
phi_33
double phi_33(double csi_x, double csi_y)
Definition: InterpSpline.c:75
obsEstimateBilin
void obsEstimateBilin(double **obsV, double *obsE, double *parV, double deltX, double deltY, int xNm, int yNm, double xMi, double yMi, int obsN)
Definition: InterpSpline.c:595
nCorrectGrad
void nCorrectGrad(double **N, double lambda, int xNum, int yNum, double deltaX, double deltaY)
Definition: InterpSpline.c:432
node_x
void node_x(double x, int *i_x, double *csi_x, double xMin, double deltaX)
Definition: InterpSpline.c:30
nCorrectLapl
void nCorrectLapl(double **N, double lambda, int xNum, int yNum, double deltaX, double deltaY)
Definition: InterpSpline.c:211
phi_3
double phi_3(double csi)
Definition: InterpSpline.c:63
dataInterpolateBilin
double dataInterpolateBilin(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:644
obsEstimateBicubic
void obsEstimateBicubic(double **obsV, double *obsE, double *parV, double deltX, double deltY, int xNm, int yNm, double xMi, double yMi, int obsN)
Definition: InterpSpline.c:467
N
#define N
Definition: e_intersect.c:917
x
#define x
phi_4
double phi_4(double csi)
Definition: InterpSpline.c:69
phi_34
double phi_34(double csi_x, double csi_y)
Definition: InterpSpline.c:81
normalDefBilin
void normalDefBilin(double **N, double *TN, double *Q, double **obsVect, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, int obsNum, int parNum, int BW)
Definition: InterpSpline.c:313
node_y
void node_y(double y, int *i_y, double *csi_y, double yMin, double deltaY)
Definition: InterpSpline.c:42
phi_44
double phi_44(double csi_x, double csi_y)
Definition: InterpSpline.c:93
normalDefBicubic
void normalDefBicubic(double **N, double *TN, double *Q, double **obsVect, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, int obsNum, int parNum, int BW)
Definition: InterpSpline.c:108
dataInterpolateBicubic
double dataInterpolateBicubic(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:533
phi_43
double phi_43(double csi_x, double csi_y)
Definition: InterpSpline.c:87
order
int order(int i_x, int i_y, int yNum)
Definition: InterpSpline.c:54
phi
double phi(double csi_x, double csi_y)
Definition: InterpSpline.c:99