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