GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
ressegm2d.c
Go to the documentation of this file.
1/*-
2 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
3 * University of Illinois
4 * US Army Construction Engineering Research Lab
5 * Copyright 1993, H. Mitasova (University of Illinois),
6 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
7 *
8 * modified by McCauley in August 1995
9 * modified by Mitasova in August 1995
10 *
11 * bug fixes by Jaro Hofierka in February 1999:
12 * line: 175,348 (*dnorm)
13 * 177,350 (points[m1].sm)
14 * 457,461 (})
15 *
16 * modified by Mitasova November 1999 (option for dnorm ind. tension)
17 *
18 */
19
20#include <stdio.h>
21#include <stdlib.h>
22#include <math.h>
23
24#include <grass/gis.h>
25#include <grass/raster.h>
26#include <grass/interpf.h>
27#include <grass/gmath.h>
28
29static int input_data(struct interp_params *, int, int, struct fcell_triple *,
30 int, int, int, int, double, double, double);
31static int write_zeros(struct interp_params *, struct quaddata *, off_t);
32
34 struct interp_params *params, struct BM *bitmask, /* bitmask */
35 double zmin, double zmax, /* min and max input z-values */
36 double *zminac, double *zmaxac, /* min and max interp. z-values */
37 double *gmin, double *gmax, /* min and max inperp. slope val. */
38 double *c1min, double *c1max, double *c2min,
39 double *c2max, /* min and max interp. curv. val. */
40 double *ertot, /* total interplating func. error */
41 off_t offset1, /* offset for temp file writing */
42 double *dnorm, int overlap, int inp_rows, int inp_cols, int fdsmooth,
43 int fdinp, double ns_res, double ew_res, double inp_ns_res,
44 double inp_ew_res, int dtens)
45{
46
47 int i, j, k, l, m, m1, i1; /* loop coounters */
48 int cursegm = 0;
49 int new_comp = 0;
50 int n_rows, n_cols /*, inp_r, inp_c */;
51 double x_or, y_or, xm, ym;
52 static int first = 1, new_first = 1;
53 double **matrix = NULL, **new_matrix = NULL, *b = NULL;
54 int *indx = NULL, *new_indx = NULL;
55 static struct fcell_triple *in_points = NULL; /* input points */
56 int out_check_rows, out_check_cols; /* total output rows/cols */
57 int first_row, last_row; /* first and last input row of segment */
58 int first_col, last_col; /* first and last input col of segment */
59 int num, prev;
60 int div; /* number of divides */
61 int rem_out_row, rem_out_col; /* output rows/cols remainders */
62 int inp_seg_r, inp_seg_c, /* # of input rows/cols in segment */
63 out_seg_r, out_seg_c; /* # of output rows/cols in segment */
64 int ngstc, nszc /* first and last output col of the
65 * segment */
66 ,
67 ngstr, nszr; /* first and last output row of the
68 * segment */
69 int index; /* index for input data */
70 int c, r;
71 int overlap1;
72 int p_size;
73 struct quaddata *data;
74 double xmax, xmin, ymax, ymin;
75 int totsegm; /* total number of segments */
76 int total_points = 0;
77 struct triple triple = {0.0, 0.0, 0.0, 0.0}; /* contains garbage */
78
79 xmin = params->x_orig;
80 ymin = params->y_orig;
81 xmax = xmin + ew_res * params->nsizc;
82 ymax = ymin + ns_res * params->nsizr;
83 prev = inp_rows * inp_cols;
84 if (prev <= params->kmax)
85 div = 1; /* no segmentation */
86
87 else { /* find the number of divides */
88 for (i = 2;; i++) {
89 c = inp_cols / i;
90 r = inp_rows / i;
91 num = c * r;
92 if (num < params->kmin) {
93 if (((params->kmin - num) > (prev + 1 - params->kmax)) &&
94 (prev + 1 < params->KMAX2)) {
95 div = i - 1;
96 break;
97 }
98 else {
99 div = i;
100 break;
101 }
102 }
103 if ((num > params->kmin) && (num + 1 < params->kmax)) {
104 div = i;
105 break;
106 }
107 prev = num;
108 }
109 }
110 out_seg_r = params->nsizr / div; /* output rows per segment */
111 out_seg_c = params->nsizc / div; /* output cols per segment */
112 inp_seg_r = inp_rows / div; /* input rows per segment */
113 inp_seg_c = inp_cols / div; /* input rows per segment */
114 rem_out_col = params->nsizc % div;
115 rem_out_row = params->nsizr % div;
116 overlap1 = min1(overlap, inp_seg_c - 1);
117 overlap1 = min1(overlap1, inp_seg_r - 1);
118 out_check_rows = 0;
119 out_check_cols = 0;
120
121 if (div == 1) {
122 p_size = inp_seg_c * inp_seg_r;
123 }
124 else {
125 p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r);
126 }
127 if (!in_points) {
128 if (!(in_points = (struct fcell_triple *)G_malloc(
129 sizeof(struct fcell_triple) * p_size * div))) {
130 fprintf(stderr, "Cannot allocate memory for in_points\n");
131 return -1;
132 }
133 }
134
135 *dnorm =
136 sqrt(((xmax - xmin) * (ymax - ymin) * p_size) / (inp_rows * inp_cols));
137
138 if (dtens) {
139 params->fi = params->fi * (*dnorm) / 1000.;
140 fprintf(stderr, "dnorm = %f, rescaled tension = %f\n", *dnorm,
141 params->fi);
142 }
143
144 if (div == 1) { /* no segmentation */
145 totsegm = 1;
146 cursegm = 1;
147
148 input_data(params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows,
149 inp_cols, zmin, inp_ns_res, inp_ew_res);
150
151 x_or = 0.;
152 y_or = 0.;
153 xm = params->nsizc * ew_res;
154 ym = params->nsizr * ns_res;
155
156 data = (struct quaddata *)quad_data_new(
157 x_or, y_or, xm, ym, params->nsizr, params->nsizc, 0, params->KMAX2);
158 m1 = 0;
159 for (k = 1; k <= p_size; k++) {
160 if (!Rast_is_f_null_value(&(in_points[k - 1].z))) {
161 data->points[m1].x = in_points[k - 1].x / (*dnorm);
162 data->points[m1].y = in_points[k - 1].y / (*dnorm);
163 /* data->points[m1].z = (double) (in_points[k - 1].z) /
164 * (*dnorm); */
165 data->points[m1].z = (double)(in_points[k - 1].z);
166 data->points[m1].sm = in_points[k - 1].smooth;
167 m1++;
168 }
169 }
170 data->n_points = m1;
171 total_points = m1;
172 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
173 fprintf(stderr, "Cannot allocate memory for indx\n");
174 return -1;
175 }
176 if (!(matrix = G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
177 fprintf(stderr, "Cannot allocate memory for matrix\n");
178 return -1;
179 }
180 if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
181 fprintf(stderr, "Cannot allocate memory for b\n");
182 return -1;
183 }
184
185 if (params->matrix_create(params, data->points, m1, matrix, indx) < 0)
186 return -1;
187 for (i = 0; i < m1; i++) {
188 b[i + 1] = data->points[i].z;
189 }
190 b[0] = 0.;
191 G_lubksb(matrix, m1 + 1, indx, b);
192
193 params->check_points(params, data, b, ertot, zmin, *dnorm, triple);
194
195 if (params->grid_calc(params, data, bitmask, zmin, zmax, zminac, zmaxac,
196 gmin, gmax, c1min, c1max, c2min, c2max, ertot, b,
197 offset1, *dnorm) < 0) {
198 fprintf(stderr, "interpolation failed\n");
199 return -1;
200 }
201 else {
202 if (totsegm != 0) {
203 G_percent(cursegm, totsegm, 1);
204 }
205 /*
206 * if (b) G_free_vector(b); if (matrix) G_free_matrix(matrix); if
207 * (indx) G_free_ivector(indx);
208 */
209 fprintf(stderr, "dnorm in ressegm after grid before out= %f \n",
210 *dnorm);
211 return total_points;
212 }
213 }
214
215 out_seg_r = params->nsizr / div; /* output rows per segment */
216 out_seg_c = params->nsizc / div; /* output cols per segment */
217 inp_seg_r = inp_rows / div; /* input rows per segment */
218 inp_seg_c = inp_cols / div; /* input rows per segment */
219 rem_out_col = params->nsizc % div;
220 rem_out_row = params->nsizr % div;
221 overlap1 = min1(overlap, inp_seg_c - 1);
222 overlap1 = min1(overlap1, inp_seg_r - 1);
223 out_check_rows = 0;
224 out_check_cols = 0;
225
226 totsegm = div * div;
227
228 /* set up a segment */
229 for (i = 1; i <= div; i++) { /* input and output rows */
230 if (i <= div - rem_out_row)
231 n_rows = out_seg_r;
232 else
233 n_rows = out_seg_r + 1;
234 /* inp_r = inp_seg_r; */
235 out_check_cols = 0;
236 ngstr = out_check_rows + 1; /* first output row of the segment */
237 nszr = ngstr + n_rows - 1; /* last output row of the segment */
238 y_or = (ngstr - 1) * ns_res; /* y origin of the segment */
239 /*
240 * Calculating input starting and ending rows and columns of this
241 * segment
242 */
243 first_row = (int)(y_or / inp_ns_res) + 1;
244 if (first_row > overlap1) {
245 first_row -= overlap1; /* middle */
246 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
247 if (last_row > inp_rows) {
248 first_row -= (last_row - inp_rows); /* bottom */
249 last_row = inp_rows;
250 }
251 }
252 else {
253 first_row = 1; /* top */
254 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
255 }
256 if ((last_row > inp_rows) || (first_row < 1)) {
257 fprintf(stderr, "Row overlap too large!\n");
258 return -1;
259 }
260 input_data(params, first_row, last_row, in_points, fdsmooth, fdinp,
261 inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res);
262
263 for (j = 1; j <= div; j++) { /* input and output cols */
264 if (j <= div - rem_out_col)
265 n_cols = out_seg_c;
266 else
267 n_cols = out_seg_c + 1;
268 /* inp_c = inp_seg_c; */
269
270 ngstc = out_check_cols + 1; /* first output col of the segment */
271 nszc = ngstc + n_cols - 1; /* last output col of the segment */
272 x_or = (ngstc - 1) * ew_res; /* x origin of the segment */
273
274 first_col = (int)(x_or / inp_ew_res) + 1;
275 if (first_col > overlap1) {
276 first_col -= overlap1; /* middle */
277 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
278 if (last_col > inp_cols) {
279 first_col -= (last_col - inp_cols); /* right */
280 last_col = inp_cols;
281 }
282 }
283 else {
284 first_col = 1; /* left */
285 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
286 }
287 if ((last_col > inp_cols) || (first_col < 1)) {
288 fprintf(stderr, "Column overlap too large!\n");
289 return -1;
290 }
291 m = 0;
292 /* Getting points for interpolation (translated) */
293
294 xm = nszc * ew_res;
295 ym = nszr * ns_res;
296 data = (struct quaddata *)quad_data_new(
297 x_or, y_or, xm, ym, nszr - ngstr + 1, nszc - ngstc + 1, 0,
298 params->KMAX2);
299 new_comp = 0;
300
301 for (k = 0; k <= last_row - first_row; k++) {
302 for (l = first_col - 1; l < last_col; l++) {
303 index = k * inp_cols + l;
304 if (!Rast_is_f_null_value(&(in_points[index].z))) {
305 /* if the point is inside the segment (not overlapping)
306 */
307 if ((in_points[index].x - x_or >= 0) &&
308 (in_points[index].y - y_or >= 0) &&
309 ((nszc - 1) * ew_res - in_points[index].x >= 0) &&
310 ((nszr - 1) * ns_res - in_points[index].y >= 0))
311 total_points += 1;
312 data->points[m].x =
313 (in_points[index].x - x_or) / (*dnorm);
314 data->points[m].y =
315 (in_points[index].y - y_or) / (*dnorm);
316 /* data->points[m].z = (double)
317 * (in_points[index].z) / (*dnorm); */
318 data->points[m].z = (double)(in_points[index].z);
319 data->points[m].sm = in_points[index].smooth;
320 m++;
321 }
322 else
323 new_comp = 1;
324
325 /* fprintf(stderr,"%f,%f,%f
326 zmin=%f\n",in_points[index].x,in_points[index].y,in_points[index].z,zmin);
327 */
328 }
329 }
330 /* fprintf (stdout,"m,index:%di,%d\n",m,index); */
331 if (m <= params->KMAX2)
332 data->n_points = m;
333 else
334 data->n_points = params->KMAX2;
335 out_check_cols += n_cols;
336 cursegm = (i - 1) * div + j - 1;
337
338 /* show before to catch 0% */
339 if (totsegm != 0) {
340 G_percent(cursegm, totsegm, 1);
341 }
342 if (m == 0) {
343 /*
344 * fprintf(stderr,"Warning: segment with zero points
345 * encountered, insrease overlap\n");
346 */
347 write_zeros(params, data, offset1);
348 }
349 else {
350 if (new_comp) {
351 if (new_first) {
352 new_first = 0;
353 if (!b) {
354 if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
355 fprintf(stderr,
356 "Cannot allocate memory for b\n");
357 return -1;
358 }
359 }
360 if (!(new_indx = G_alloc_ivector(params->KMAX2 + 1))) {
361 fprintf(stderr,
362 "Cannot allocate memory for new_indx\n");
363 return -1;
364 }
365 if (!(new_matrix = G_alloc_matrix(params->KMAX2 + 1,
366 params->KMAX2 + 1))) {
367 fprintf(stderr,
368 "Cannot allocate memory for new_matrix\n");
369 return -1;
370 }
371 } /*new_first */
372 if (params->matrix_create(params, data->points,
373 data->n_points, new_matrix,
374 new_indx) < 0)
375 return -1;
376
377 for (i1 = 0; i1 < m; i1++) {
378 b[i1 + 1] = data->points[i1].z;
379 }
380 b[0] = 0.;
381 G_lubksb(new_matrix, data->n_points + 1, new_indx, b);
382
383 params->check_points(params, data, b, ertot, zmin, *dnorm,
384 triple);
385
386 if (params->grid_calc(params, data, bitmask, zmin, zmax,
387 zminac, zmaxac, gmin, gmax, c1min,
388 c1max, c2min, c2max, ertot, b,
389 offset1, *dnorm) < 0) {
390
391 fprintf(stderr, "interpolate() failed\n");
392 return -1;
393 }
394 } /*new_comp */
395 else {
396 if (first) {
397 first = 0;
398 if (!b) {
399 if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
400 fprintf(stderr,
401 "Cannot allocate memory for b\n");
402 return -1;
403 }
404 }
405 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
406 fprintf(stderr,
407 "Cannot allocate memory for indx\n");
408 return -1;
409 }
410 if (!(matrix = G_alloc_matrix(params->KMAX2 + 1,
411 params->KMAX2 + 1))) {
412 fprintf(stderr,
413 "Cannot allocate memory for matrix\n");
414 return -1;
415 }
416 } /* first */
417 if (params->matrix_create(params, data->points,
418 data->n_points, matrix, indx) < 0)
419 return -1;
420 /* } here it was bug */
421 for (i1 = 0; i1 < m; i1++)
422 b[i1 + 1] = data->points[i1].z;
423 b[0] = 0.;
424 G_lubksb(matrix, data->n_points + 1, indx, b);
425
426 params->check_points(params, data, b, ertot, zmin, *dnorm,
427 triple);
428
429 if (params->grid_calc(params, data, bitmask, zmin, zmax,
430 zminac, zmaxac, gmin, gmax, c1min,
431 c1max, c2min, c2max, ertot, b,
432 offset1, *dnorm) < 0) {
433
434 fprintf(stderr, "interpolate() failed\n");
435 return -1;
436 }
437 }
438 }
439 if (data) {
440 G_free(data->points);
441 G_free(data);
442 }
443 /*
444 * cursegm++;
445 */
446 }
447 out_check_rows += n_rows;
448 }
449
450 /* run one last time after the loop is done to catch 100% */
451 if (totsegm != 0)
452 G_percent(1, 1,
453 1); /* cursegm doesn't get to totsegm so we force 100% */
454
455 /*
456 * if (b) G_free_vector(b); if (indx) G_free_ivector(indx); if (matrix)
457 * G_free_matrix(matrix);
458 */
459 fprintf(stderr, "dnorm in ressegm after grid before out2= %f \n", *dnorm);
460 return total_points;
461}
462
463/* input of data for interpolation and smoothing parameters */
464
465static int input_data(struct interp_params *params, int first_row, int last_row,
466 struct fcell_triple *points, int fdsmooth, int fdinp,
467 int inp_rows, int inp_cols, double zmin,
468 double inp_ns_res, double inp_ew_res)
469{
470 double x, y, sm; /* input data and smoothing */
471 int m1, m2; /* loop counters */
472 static FCELL *cellinp = NULL; /* cell buffer for input data */
473 static FCELL *cellsmooth = NULL; /* cell buffer for smoothing */
474
475 if (!cellinp)
476 cellinp = Rast_allocate_f_buf();
477 if (!cellsmooth)
478 cellsmooth = Rast_allocate_f_buf();
479
480 for (m1 = 0; m1 <= last_row - first_row; m1++) {
481 Rast_get_f_row(fdinp, cellinp, inp_rows - m1 - first_row);
482 if (fdsmooth >= 0)
483 Rast_get_f_row(fdsmooth, cellsmooth, inp_rows - m1 - first_row);
484
485 y = params->y_orig + (m1 + first_row - 1 + 0.5) * inp_ns_res;
486 for (m2 = 0; m2 < inp_cols; m2++) {
487 x = params->x_orig + (m2 + 0.5) * inp_ew_res;
488 /*
489 * z = cellinp[m2]*params->zmult;
490 */
491 if (fdsmooth >= 0)
492 sm = (double)cellsmooth[m2];
493 else
494 sm = 0.01;
495
496 points[m1 * inp_cols + m2].x = x - params->x_orig;
497 points[m1 * inp_cols + m2].y = y - params->y_orig;
498 if (!Rast_is_f_null_value(cellinp + m2)) {
499 points[m1 * inp_cols + m2].z =
500 cellinp[m2] * params->zmult - zmin;
501 }
502 else {
503 Rast_set_f_null_value(&(points[m1 * inp_cols + m2].z), 1);
504 }
505
506 /* fprintf (stdout,"sm: %f\n",sm); */
507
508 points[m1 * inp_cols + m2].smooth = sm;
509 }
510 }
511 return 1;
512}
513
514static int write_zeros(struct interp_params *params,
515 struct quaddata *data, /* given segment */
516 off_t offset1 /* offset for temp file writing */
517)
518{
519
520 /*
521 * C C INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
522 * c
523 */
524 double x_or = data->x_orig;
525 double y_or = data->y_orig;
526 int n_rows = data->n_rows;
527 int n_cols = data->n_cols;
528 int cond1, cond2;
529 int k, l;
530 int ngstc, nszc, ngstr, nszr;
531 off_t offset, offset2;
532 double ns_res, ew_res;
533
534 ns_res = (((struct quaddata *)(data))->ymax -
535 ((struct quaddata *)(data))->y_orig) /
536 data->n_rows;
537 ew_res = (((struct quaddata *)(data))->xmax -
538 ((struct quaddata *)(data))->x_orig) /
539 data->n_cols;
540
541 cond2 = ((params->adxx != NULL) || (params->adyy != NULL) ||
542 (params->adxy != NULL));
543 cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);
544
545 ngstc = (int)(x_or / ew_res + 0.5) + 1;
546 nszc = ngstc + n_cols - 1;
547 ngstr = (int)(y_or / ns_res + 0.5) + 1;
548 nszr = ngstr + n_rows - 1;
549
550 for (k = ngstr; k <= nszr; k++) {
551 offset = offset1 * (k - 1); /* rows offset */
552 for (l = ngstc; l <= nszc; l++) {
553 /*
554 * params->az[l] = 0.;
555 */
556 Rast_set_d_null_value(params->az + l, 1);
557 if (cond1) {
558 /*
559 * params->adx[l] = (FCELL)0.; params->ady[l] = (FCELL)0.;
560 */
561 Rast_set_d_null_value(params->adx + l, 1);
562 Rast_set_d_null_value(params->ady + l, 1);
563 if (cond2) {
564 Rast_set_d_null_value(params->adxx + l, 1);
565 Rast_set_d_null_value(params->adyy + l, 1);
566 Rast_set_d_null_value(params->adxy + l, 1);
567 /*
568 * params->adxx[l] = (FCELL)0.; params->adyy[l] = (FCELL)0.;
569 * params->adxy[l] = (FCELL)0.;
570 */
571 }
572 }
573 }
574 offset2 = (offset + ngstc - 1) * sizeof(FCELL);
575 if (params->wr_temp(params, ngstc, nszc, offset2) < 0)
576 return -1;
577 }
578 return 1;
579}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#define NULL
Definition ccmath.h:32
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:57
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:39
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
Definition dataquad.c:59
double b
double l
double r
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition ialloc.c:38
int min1(int, int)
Definition minmax.c:17
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition lu.c:103
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
int IL_resample_interp_segments_2d(struct interp_params *params, struct BM *bitmask, double zmin, double zmax, double *zminac, double *zmaxac, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, double *ertot, off_t offset1, double *dnorm, int overlap, int inp_rows, int inp_cols, int fdsmooth, int fdinp, double ns_res, double ew_res, double inp_ns_res, double inp_ew_res, int dtens)
Definition ressegm2d.c:33
double x
Definition interpf.h:18
double smooth
Definition interpf.h:21
double y
Definition interpf.h:19
check_points_fn * check_points
Definition interpf.h:120
double zmult
Definition interpf.h:67
DCELL * az
Definition interpf.h:85
grid_calc_fn * grid_calc
Definition interpf.h:116
double fi
Definition interpf.h:88
double x_orig
Definition interpf.h:101
DCELL * adxy
Definition interpf.h:85
DCELL * adyy
Definition interpf.h:85
DCELL * adx
Definition interpf.h:85
double y_orig
Definition interpf.h:101
DCELL * ady
Definition interpf.h:85
wr_temp_fn * wr_temp
Definition interpf.h:128
DCELL * adxx
Definition interpf.h:85
matrix_create_fn * matrix_create
Definition interpf.h:118
double ymax
Definition dataquad.h:49
double y_orig
Definition dataquad.h:47
double x_orig
Definition dataquad.h:46
struct triple * points
Definition dataquad.h:53
int n_points
Definition dataquad.h:52
double xmax
Definition dataquad.h:48
int n_cols
Definition dataquad.h:51
int n_rows
Definition dataquad.h:50
double z
Definition dataquad.h:41
double sm
Definition dataquad.h:42
double x
Definition dataquad.h:39
double y
Definition dataquad.h:40
#define x