GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993 00004 * University of Illinois 00005 * US Army Construction Engineering Research Lab 00006 * Copyright 1993, H. Mitasova (University of Illinois), 00007 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL) 00008 * 00009 * modified by McCauley in August 1995 00010 * modified by Mitasova in August 1995 00011 * 00012 * bug fixes by Jaro Hofierka in February 1999: 00013 * line: 175,348 (*dnorm) 00014 * 177,350 (points[m1].sm) 00015 * 457,461 (}) 00016 * 00017 * modified by Mitasova November 1999 (option for dnorm ind. tension) 00018 * 00019 */ 00020 00021 #include <stdio.h> 00022 #include <stdlib.h> 00023 #include <math.h> 00024 #include <grass/gis.h> 00025 #include <grass/interpf.h> 00026 #include <grass/gmath.h> 00027 00028 static int input_data(struct interp_params *, 00029 int, int, struct fcell_triple *, int, int, int, int, 00030 double, double, double); 00031 static int write_zeros(struct interp_params *, struct quaddata *, int); 00032 00033 int IL_resample_interp_segments_2d(struct interp_params *params, struct BM *bitmask, /* bitmask */ 00034 double zmin, double zmax, /* min and max input z-values */ 00035 double *zminac, double *zmaxac, /* min and max interp. z-values */ 00036 double *gmin, double *gmax, /* min and max inperp. slope val. */ 00037 double *c1min, double *c1max, double *c2min, double *c2max, /* min and max interp. curv. val. */ 00038 double *ertot, /* total interplating func. error */ 00039 int offset1, /* offset for temp file writing */ 00040 double *dnorm, 00041 int overlap, 00042 int inp_rows, 00043 int inp_cols, 00044 int fdsmooth, 00045 int fdinp, 00046 double ns_res, 00047 double ew_res, 00048 double inp_ns_res, 00049 double inp_ew_res, int dtens) 00050 { 00051 00052 int i, j, k, l, m, m1, i1; /* loop coounters */ 00053 int cursegm = 0; 00054 int new_comp = 0; 00055 int n_rows, n_cols, inp_r, inp_c; 00056 double x_or, y_or, xm, ym; 00057 static int first = 1, new_first = 1; 00058 double **matrix = NULL, **new_matrix = NULL, *b = NULL; 00059 int *indx = NULL, *new_indx = NULL; 00060 static struct fcell_triple *in_points = NULL; /* input points */ 00061 int inp_check_rows, inp_check_cols, /* total input rows/cols */ 00062 out_check_rows, out_check_cols; /* total output rows/cols */ 00063 int first_row, last_row; /* first and last input row of segment */ 00064 int first_col, last_col; /* first and last input col of segment */ 00065 int num, prev; 00066 int div; /* number of divides */ 00067 int rem_out_row, rem_out_col; /* output rows/cols remainders */ 00068 int inp_seg_r, inp_seg_c, /* # of input rows/cols in segment */ 00069 out_seg_r, out_seg_c; /* # of output rows/cols in segment */ 00070 int ngstc, nszc /* first and last output col of the 00071 * segment */ 00072 , ngstr, nszr; /* first and last output row of the 00073 * segment */ 00074 int index; /* index for input data */ 00075 int c, r; 00076 int overlap1; 00077 int p_size; 00078 struct quaddata *data; 00079 double xmax, xmin, ymax, ymin; 00080 int totsegm; /* total number of segments */ 00081 int total_points = 0; 00082 00083 00084 xmin = params->x_orig; 00085 ymin = params->y_orig; 00086 xmax = xmin + ew_res * params->nsizc; 00087 ymax = ymin + ns_res * params->nsizr; 00088 prev = inp_rows * inp_cols; 00089 if (prev <= params->kmax) 00090 div = 1; /* no segmentation */ 00091 00092 else { /* find the number of divides */ 00093 for (i = 2;; i++) { 00094 c = inp_cols / i; 00095 r = inp_rows / i; 00096 num = c * r; 00097 if (num < params->kmin) { 00098 if (((params->kmin - num) > (prev + 1 - params->kmax)) && 00099 (prev + 1 < params->KMAX2)) { 00100 div = i - 1; 00101 break; 00102 } 00103 else { 00104 div = i; 00105 break; 00106 } 00107 } 00108 if ((num > params->kmin) && (num + 1 < params->kmax)) { 00109 div = i; 00110 break; 00111 } 00112 prev = num; 00113 } 00114 } 00115 out_seg_r = params->nsizr / div; /* output rows per segment */ 00116 out_seg_c = params->nsizc / div; /* output cols per segment */ 00117 inp_seg_r = inp_rows / div; /* input rows per segment */ 00118 inp_seg_c = inp_cols / div; /* input rows per segment */ 00119 rem_out_col = params->nsizc % div; 00120 rem_out_row = params->nsizr % div; 00121 overlap1 = min1(overlap, inp_seg_c - 1); 00122 overlap1 = min1(overlap1, inp_seg_r - 1); 00123 out_check_rows = 0; 00124 out_check_cols = 0; 00125 inp_check_rows = 0; 00126 inp_check_cols = 0; 00127 00128 if (div == 1) { 00129 p_size = inp_seg_c * inp_seg_r; 00130 } 00131 else { 00132 p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r); 00133 } 00134 if (!in_points) { 00135 if (! 00136 (in_points = 00137 (struct fcell_triple *)G_malloc(sizeof(struct fcell_triple) * 00138 p_size * div))) { 00139 fprintf(stderr, "Cannot allocate memory for in_points\n"); 00140 return -1; 00141 } 00142 } 00143 00144 *dnorm = 00145 sqrt(((xmax - xmin) * (ymax - 00146 ymin) * p_size) / (inp_rows * inp_cols)); 00147 00148 if (dtens) { 00149 params->fi = params->fi * (*dnorm) / 1000.; 00150 fprintf(stderr, "dnorm = %f, rescaled tension = %f\n", *dnorm, 00151 params->fi); 00152 } 00153 00154 if (div == 1) { /* no segmentation */ 00155 totsegm = 1; 00156 cursegm = 1; 00157 00158 input_data(params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows, 00159 inp_cols, zmin, inp_ns_res, inp_ew_res); 00160 00161 x_or = 0.; 00162 y_or = 0.; 00163 xm = params->nsizc * ew_res; 00164 ym = params->nsizr * ns_res; 00165 00166 data = (struct quaddata *)quad_data_new(x_or, y_or, xm, ym, 00167 params->nsizr, params->nsizc, 00168 0, params->KMAX2); 00169 m1 = 0; 00170 for (k = 1; k <= p_size; k++) { 00171 if (!G_is_f_null_value(&(in_points[k - 1].z))) { 00172 data->points[m1].x = in_points[k - 1].x / (*dnorm); 00173 data->points[m1].y = in_points[k - 1].y / (*dnorm); 00174 /* data->points[m1].z = (double) (in_points[k - 1].z) / (*dnorm); */ 00175 data->points[m1].z = (double)(in_points[k - 1].z); 00176 data->points[m1].sm = in_points[k - 1].smooth; 00177 m1++; 00178 } 00179 } 00180 data->n_points = m1; 00181 total_points = m1; 00182 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) { 00183 fprintf(stderr, "Cannot allocate memory for indx\n"); 00184 return -1; 00185 } 00186 if (!(matrix = G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) { 00187 fprintf(stderr, "Cannot allocate memory for matrix\n"); 00188 return -1; 00189 } 00190 if (!(b = G_alloc_vector(params->KMAX2 + 2))) { 00191 fprintf(stderr, "Cannot allocate memory for b\n"); 00192 return -1; 00193 } 00194 00195 if (params->matrix_create(params, data->points, m1, matrix, indx) < 0) 00196 return -1; 00197 for (i = 0; i < m1; i++) { 00198 b[i + 1] = data->points[i].z; 00199 } 00200 b[0] = 0.; 00201 G_lubksb(matrix, m1 + 1, indx, b); 00202 00203 params->check_points(params, data, b, ertot, zmin, *dnorm); 00204 00205 if (params->grid_calc(params, data, bitmask, 00206 zmin, zmax, zminac, zmaxac, gmin, gmax, 00207 c1min, c1max, c2min, c2max, ertot, b, offset1, 00208 *dnorm) < 0) { 00209 fprintf(stderr, "interpolation failed\n"); 00210 return -1; 00211 } 00212 else { 00213 if (totsegm != 0) { 00214 G_percent(cursegm, totsegm, 1); 00215 } 00216 /* 00217 * if (b) G_free_vector(b); if (matrix) G_free_matrix(matrix); if 00218 * (indx) G_free_ivector(indx); 00219 */ 00220 fprintf(stderr, "dnorm in ressegm after grid before out= %f \n", 00221 *dnorm); 00222 return total_points; 00223 } 00224 } 00225 00226 out_seg_r = params->nsizr / div; /* output rows per segment */ 00227 out_seg_c = params->nsizc / div; /* output cols per segment */ 00228 inp_seg_r = inp_rows / div; /* input rows per segment */ 00229 inp_seg_c = inp_cols / div; /* input rows per segment */ 00230 rem_out_col = params->nsizc % div; 00231 rem_out_row = params->nsizr % div; 00232 overlap1 = min1(overlap, inp_seg_c - 1); 00233 overlap1 = min1(overlap1, inp_seg_r - 1); 00234 out_check_rows = 0; 00235 out_check_cols = 0; 00236 inp_check_rows = 0; 00237 inp_check_cols = 0; 00238 00239 totsegm = div * div; 00240 00241 /* set up a segment */ 00242 for (i = 1; i <= div; i++) { /* input and output rows */ 00243 if (i <= div - rem_out_row) 00244 n_rows = out_seg_r; 00245 else 00246 n_rows = out_seg_r + 1; 00247 inp_r = inp_seg_r; 00248 out_check_cols = 0; 00249 inp_check_cols = 0; 00250 ngstr = out_check_rows + 1; /* first output row of the segment */ 00251 nszr = ngstr + n_rows - 1; /* last output row of the segment */ 00252 y_or = (ngstr - 1) * ns_res; /* y origin of the segment */ 00253 /* 00254 * Calculating input starting and ending rows and columns of this 00255 * segment 00256 */ 00257 first_row = (int)(y_or / inp_ns_res) + 1; 00258 if (first_row > overlap1) { 00259 first_row -= overlap1; /* middle */ 00260 last_row = first_row + inp_seg_r + overlap1 * 2 - 1; 00261 if (last_row > inp_rows) { 00262 first_row -= (last_row - inp_rows); /* bottom */ 00263 last_row = inp_rows; 00264 } 00265 } 00266 else { 00267 first_row = 1; /* top */ 00268 last_row = first_row + inp_seg_r + overlap1 * 2 - 1; 00269 } 00270 if ((last_row > inp_rows) || (first_row < 1)) { 00271 fprintf(stderr, "Row overlap too large!\n"); 00272 return -1; 00273 } 00274 input_data(params, first_row, last_row, in_points, fdsmooth, fdinp, 00275 inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res); 00276 00277 for (j = 1; j <= div; j++) { /* input and output cols */ 00278 if (j <= div - rem_out_col) 00279 n_cols = out_seg_c; 00280 else 00281 n_cols = out_seg_c + 1; 00282 inp_c = inp_seg_c; 00283 00284 ngstc = out_check_cols + 1; /* first output col of the segment */ 00285 nszc = ngstc + n_cols - 1; /* last output col of the segment */ 00286 x_or = (ngstc - 1) * ew_res; /* x origin of the segment */ 00287 00288 first_col = (int)(x_or / inp_ew_res) + 1; 00289 if (first_col > overlap1) { 00290 first_col -= overlap1; /* middle */ 00291 last_col = first_col + inp_seg_c + overlap1 * 2 - 1; 00292 if (last_col > inp_cols) { 00293 first_col -= (last_col - inp_cols); /* right */ 00294 last_col = inp_cols; 00295 } 00296 } 00297 else { 00298 first_col = 1; /* left */ 00299 last_col = first_col + inp_seg_c + overlap1 * 2 - 1; 00300 } 00301 if ((last_col > inp_cols) || (first_col < 1)) { 00302 fprintf(stderr, "Column overlap too large!\n"); 00303 return -1; 00304 } 00305 m = 0; 00306 /* Getting points for interpolation (translated) */ 00307 00308 xm = nszc * ew_res; 00309 ym = nszr * ns_res; 00310 data = (struct quaddata *)quad_data_new(x_or, y_or, xm, ym, 00311 nszr - ngstr + 1, 00312 nszc - ngstc + 1, 0, 00313 params->KMAX2); 00314 new_comp = 0; 00315 00316 for (k = 0; k <= last_row - first_row; k++) { 00317 for (l = first_col - 1; l < last_col; l++) { 00318 index = k * inp_cols + l; 00319 if (!G_is_f_null_value(&(in_points[index].z))) { 00320 /* if the point is inside the segment (not overlapping) */ 00321 if ((in_points[index].x - x_or >= 0) && 00322 (in_points[index].y - y_or >= 0) && 00323 ((nszc - 1) * ew_res - in_points[index].x >= 0) && 00324 ((nszr - 1) * ns_res - in_points[index].y >= 0)) 00325 total_points += 1; 00326 data->points[m].x = 00327 (in_points[index].x - x_or) / (*dnorm); 00328 data->points[m].y = 00329 (in_points[index].y - y_or) / (*dnorm); 00330 /* data->points[m].z = (double) (in_points[index].z) / (*dnorm); */ 00331 data->points[m].z = (double)(in_points[index].z); 00332 data->points[m].sm = in_points[index].smooth; 00333 m++; 00334 } 00335 else 00336 new_comp = 1; 00337 00338 /* fprintf(stderr,"%f,%f,%f 00339 zmin=%f\n",in_points[index].x,in_points[index].y,in_points[index].z,zmin); 00340 */ 00341 } 00342 } 00343 /* fprintf (stdout,"m,index:%di,%d\n",m,index); */ 00344 if (m <= params->KMAX2) 00345 data->n_points = m; 00346 else 00347 data->n_points = params->KMAX2; 00348 out_check_cols += n_cols; 00349 inp_check_cols += inp_c; 00350 cursegm = (i - 1) * div + j - 1; 00351 00352 /* show before to catch 0% */ 00353 if (totsegm != 0) { 00354 G_percent(cursegm, totsegm, 1); 00355 } 00356 if (m == 0) { 00357 /* 00358 * fprintf(stderr,"Warning: segment with zero points encountered, 00359 * insrease overlap\n"); 00360 */ 00361 write_zeros(params, data, offset1); 00362 } 00363 else { 00364 if (new_comp) { 00365 if (new_first) { 00366 new_first = 0; 00367 if (!b) { 00368 if (!(b = G_alloc_vector(params->KMAX2 + 2))) { 00369 fprintf(stderr, 00370 "Cannot allocate memory for b\n"); 00371 return -1; 00372 } 00373 } 00374 if (!(new_indx = G_alloc_ivector(params->KMAX2 + 1))) { 00375 fprintf(stderr, 00376 "Cannot allocate memory for new_indx\n"); 00377 return -1; 00378 } 00379 if (! 00380 (new_matrix = 00381 G_alloc_matrix(params->KMAX2 + 1, 00382 params->KMAX2 + 1))) { 00383 fprintf(stderr, 00384 "Cannot allocate memory for new_matrix\n"); 00385 return -1; 00386 } 00387 } /*new_first */ 00388 if (params-> 00389 matrix_create(params, data->points, data->n_points, 00390 new_matrix, new_indx) < 0) 00391 return -1; 00392 00393 for (i1 = 0; i1 < m; i1++) { 00394 b[i1 + 1] = data->points[i1].z; 00395 } 00396 b[0] = 0.; 00397 G_lubksb(new_matrix, data->n_points + 1, new_indx, b); 00398 00399 params->check_points(params, data, b, ertot, zmin, 00400 *dnorm); 00401 00402 if (params->grid_calc(params, data, bitmask, 00403 zmin, zmax, zminac, zmaxac, gmin, 00404 gmax, c1min, c1max, c2min, c2max, 00405 ertot, b, offset1, *dnorm) < 0) { 00406 00407 fprintf(stderr, "interpolate() failed\n"); 00408 return -1; 00409 } 00410 } /*new_comp */ 00411 else { 00412 if (first) { 00413 first = 0; 00414 if (!b) { 00415 if (!(b = G_alloc_vector(params->KMAX2 + 2))) { 00416 fprintf(stderr, 00417 "Cannot allocate memory for b\n"); 00418 return -1; 00419 } 00420 } 00421 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) { 00422 fprintf(stderr, 00423 "Cannot allocate memory for indx\n"); 00424 return -1; 00425 } 00426 if (! 00427 (matrix = 00428 G_alloc_matrix(params->KMAX2 + 1, 00429 params->KMAX2 + 1))) { 00430 fprintf(stderr, 00431 "Cannot allocate memory for matrix\n"); 00432 return -1; 00433 } 00434 } /* first */ 00435 if (params-> 00436 matrix_create(params, data->points, data->n_points, 00437 matrix, indx) < 0) 00438 return -1; 00439 /* } here it was bug */ 00440 for (i1 = 0; i1 < m; i1++) 00441 b[i1 + 1] = data->points[i1].z; 00442 b[0] = 0.; 00443 G_lubksb(matrix, data->n_points + 1, indx, b); 00444 00445 params->check_points(params, data, b, ertot, zmin, 00446 *dnorm); 00447 00448 if (params->grid_calc(params, data, bitmask, 00449 zmin, zmax, zminac, zmaxac, gmin, 00450 gmax, c1min, c1max, c2min, c2max, 00451 ertot, b, offset1, *dnorm) < 0) { 00452 00453 fprintf(stderr, "interpolate() failed\n"); 00454 return -1; 00455 } 00456 } 00457 } 00458 if (data) { 00459 G_free(data->points); 00460 G_free(data); 00461 } 00462 /* 00463 * cursegm++; 00464 */ 00465 } 00466 00467 inp_check_rows += inp_r; 00468 out_check_rows += n_rows; 00469 } 00470 00471 /* run one last time after the loop is done to catch 100% */ 00472 if (totsegm != 0) 00473 G_percent(1, 1, 1); /* cursegm doesn't get to totsegm so we force 100% */ 00474 00475 /* 00476 * if (b) G_free_vector(b); if (indx) G_free_ivector(indx); if (matrix) 00477 * G_free_matrix(matrix); 00478 */ 00479 fprintf(stderr, "dnorm in ressegm after grid before out2= %f \n", *dnorm); 00480 return total_points; 00481 } 00482 00483 /* input of data for interpolation and smoothing parameters */ 00484 00485 static int input_data(struct interp_params *params, 00486 int first_row, int last_row, 00487 struct fcell_triple *points, 00488 int fdsmooth, int fdinp, 00489 int inp_rows, int inp_cols, 00490 double zmin, double inp_ns_res, double inp_ew_res) 00491 { 00492 double x, y, sm; /* input data and smoothing */ 00493 int m1, m2; /* loop counters */ 00494 int ret_val, ret_val1; /* return values of G_get_map_row */ 00495 static FCELL *cellinp = NULL; /* cell buffer for input data */ 00496 static FCELL *cellsmooth = NULL; /* cell buffer for smoothing */ 00497 00498 00499 if (!cellinp) 00500 cellinp = G_allocate_f_raster_buf(); 00501 if (!cellsmooth) 00502 cellsmooth = G_allocate_f_raster_buf(); 00503 00504 for (m1 = 0; m1 <= last_row - first_row; m1++) { 00505 ret_val = 00506 G_get_f_raster_row(fdinp, cellinp, inp_rows - m1 - first_row); 00507 if (ret_val < 0) { 00508 fprintf(stderr, "Cannot get row %d (return value = %d)\n", m1, 00509 ret_val); 00510 return -1; 00511 } 00512 if (fdsmooth >= 0) { 00513 ret_val1 = 00514 G_get_f_raster_row(fdsmooth, cellsmooth, 00515 inp_rows - m1 - first_row); 00516 if (ret_val1 < 0) { 00517 fprintf(stderr, "Cannot get smoothing row\n"); 00518 } 00519 } 00520 y = params->y_orig + (m1 + first_row - 1 + 0.5) * inp_ns_res; 00521 for (m2 = 0; m2 < inp_cols; m2++) { 00522 x = params->x_orig + (m2 + 0.5) * inp_ew_res; 00523 /* 00524 * z = cellinp[m2]*params->zmult; 00525 */ 00526 if (fdsmooth >= 0) 00527 sm = (double)cellsmooth[m2]; 00528 else 00529 sm = 0.01; 00530 00531 points[m1 * inp_cols + m2].x = x - params->x_orig; 00532 points[m1 * inp_cols + m2].y = y - params->y_orig; 00533 if (!G_is_f_null_value(cellinp + m2)) { 00534 points[m1 * inp_cols + m2].z = 00535 cellinp[m2] * params->zmult - zmin; 00536 } 00537 else { 00538 G_set_f_null_value(&(points[m1 * inp_cols + m2].z), 1); 00539 } 00540 00541 /* fprintf (stdout,"sm: %f\n",sm); */ 00542 00543 points[m1 * inp_cols + m2].smooth = sm; 00544 } 00545 } 00546 return 1; 00547 } 00548 00549 static int write_zeros(struct interp_params *params, struct quaddata *data, /* given segment */ 00550 int offset1 /* offset for temp file writing */ 00551 ) 00552 { 00553 00554 /* 00555 * C C INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul. 00556 * c 00557 */ 00558 double x_or = data->x_orig; 00559 double y_or = data->y_orig; 00560 int n_rows = data->n_rows; 00561 int n_cols = data->n_cols; 00562 int cond1, cond2; 00563 int k, l; 00564 int ngstc, nszc, ngstr, nszr; 00565 int offset, offset2; 00566 double ns_res, ew_res; 00567 00568 ns_res = (((struct quaddata *)(data))->ymax - 00569 ((struct quaddata *)(data))->y_orig) / data->n_rows; 00570 ew_res = (((struct quaddata *)(data))->xmax - 00571 ((struct quaddata *)(data))->x_orig) / data->n_cols; 00572 00573 cond2 = ((params->adxx != NULL) || (params->adyy != NULL) || 00574 (params->adxy != NULL)); 00575 cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2); 00576 00577 ngstc = (int)(x_or / ew_res + 0.5) + 1; 00578 nszc = ngstc + n_cols - 1; 00579 ngstr = (int)(y_or / ns_res + 0.5) + 1; 00580 nszr = ngstr + n_rows - 1; 00581 00582 for (k = ngstr; k <= nszr; k++) { 00583 offset = offset1 * (k - 1); /* rows offset */ 00584 for (l = ngstc; l <= nszc; l++) { 00585 /* 00586 * params->az[l] = 0.; 00587 */ 00588 G_set_d_null_value(params->az + l, 1); 00589 if (cond1) { 00590 /* 00591 * params->adx[l] = (FCELL)0.; params->ady[l] = (FCELL)0.; 00592 */ 00593 G_set_d_null_value(params->adx + l, 1); 00594 G_set_d_null_value(params->ady + l, 1); 00595 if (cond2) { 00596 G_set_d_null_value(params->adxx + l, 1); 00597 G_set_d_null_value(params->adyy + l, 1); 00598 G_set_d_null_value(params->adxy + l, 1); 00599 /* 00600 * params->adxx[l] = (FCELL)0.; params->adyy[l] = (FCELL)0.; 00601 * params->adxy[l] = (FCELL)0.; 00602 */ 00603 } 00604 } 00605 } 00606 offset2 = (offset + ngstc - 1) * sizeof(FCELL); 00607 if (params->wr_temp(params, ngstc, nszc, offset2) < 0) 00608 return -1; 00609 } 00610 return 1; 00611 }