GRASS Programmer's Manual
6.4.2(2012)
|
00001 /* 00002 ** Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 1993 00003 ** Copyright H. Mitasova, I. Kosinovsky, D.Gerdes 00004 */ 00005 00006 #include <stdio.h> 00007 #include <stdlib.h> 00008 #include <math.h> 00009 #include <grass/gis.h> 00010 #include <grass/glocale.h> 00011 #include <grass/interpf.h> 00012 #include <grass/gmath.h> 00013 00014 static double smallest_segment(struct multtree *, int); 00015 00016 int IL_interp_segments_2d(struct interp_params *params, struct tree_info *info, /* info for the quad tree */ 00017 struct multtree *tree, /* current leaf of the quad tree */ 00018 struct BM *bitmask, /* bitmask */ 00019 double zmin, double zmax, /* min and max input z-values */ 00020 double *zminac, double *zmaxac, /* min and max interp. z-values */ 00021 double *gmin, double *gmax, /* min and max inperp. slope val. */ 00022 double *c1min, double *c1max, double *c2min, double *c2max, /* min and max interp. curv. val. */ 00023 double *ertot, /* total interplating func. error */ 00024 int totsegm, /* total number of segments */ 00025 int offset1, /* offset for temp file writing */ 00026 double dnorm) 00027 /* 00028 Recursively processes each segment in a tree by 00029 a) finding points from neighbouring segments so that the total number of 00030 points is between KMIN and KMAX2 by calling tree function MT_get_region(). 00031 b) creating and solving the system of linear equations using these points 00032 and interp() by calling matrix_create() and G_ludcmp(). 00033 c) checking the interpolating function values at points by calling 00034 check_points(). 00035 d) computing grid for this segment using points and interp() by calling 00036 grid_calc(). 00037 */ 00038 { 00039 double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1, temp2; 00040 int i, npt, nptprev, MAXENC; 00041 struct quaddata *data; 00042 static int cursegm = 0; 00043 static double *b = NULL; 00044 static int *indx = NULL; 00045 static double **matrix = NULL; 00046 double ew_res, ns_res; 00047 static int first_time = 1; 00048 static double smseg; 00049 int MINPTS; 00050 double pr; 00051 struct triple *point; 00052 struct triple skip_point; 00053 int m_skip, skip_index, j, k, segtest; 00054 double xx, yy, zz; 00055 00056 /* find the size of the smallest segment once */ 00057 if (first_time) { 00058 smseg = smallest_segment(info->root, 4); 00059 first_time = 0; 00060 } 00061 ns_res = (((struct quaddata *)(info->root->data))->ymax - 00062 ((struct quaddata *)(info->root->data))->y_orig) / 00063 params->nsizr; 00064 ew_res = 00065 (((struct quaddata *)(info->root->data))->xmax - 00066 ((struct quaddata *)(info->root->data))->x_orig) / params->nsizc; 00067 00068 if (tree == NULL) 00069 return -1; 00070 if (tree->data == NULL) 00071 return -1; 00072 if (((struct quaddata *)(tree->data))->points == NULL) { 00073 for (i = 0; i < 4; i++) { 00074 IL_interp_segments_2d(params, info, tree->leafs[i], 00075 bitmask, zmin, zmax, zminac, zmaxac, gmin, 00076 gmax, c1min, c1max, c2min, c2max, ertot, 00077 totsegm, offset1, dnorm); 00078 } 00079 return 1; 00080 } 00081 else { 00082 distx = (((struct quaddata *)(tree->data))->n_cols * ew_res) * 0.1; 00083 disty = (((struct quaddata *)(tree->data))->n_rows * ns_res) * 0.1; 00084 distxp = 0; 00085 distyp = 0; 00086 xmn = ((struct quaddata *)(tree->data))->x_orig; 00087 xmx = ((struct quaddata *)(tree->data))->xmax; 00088 ymn = ((struct quaddata *)(tree->data))->y_orig; 00089 ymx = ((struct quaddata *)(tree->data))->ymax; 00090 i = 0; 00091 MAXENC = 0; 00092 /* data is a window with zero points; some fields don't make sence in this case 00093 so they are zero (like resolution,dimentions */ 00094 /* CHANGE */ 00095 /* Calcutaing kmin for surrent segment (depends on the size) */ 00096 00097 /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/ 00098 pr = pow(2., (xmx - xmn) / smseg - 1.); 00099 MINPTS = 00100 params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2)); 00101 /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */ 00102 00103 data = 00104 (struct quaddata *)quad_data_new(xmn - distx, ymn - disty, 00105 xmx + distx, ymx + disty, 0, 0, 00106 0, params->KMAX2); 00107 npt = MT_region_data(info, info->root, data, params->KMAX2, 4); 00108 00109 while ((npt < MINPTS) || (npt > params->KMAX2)) { 00110 if (i >= 70) { 00111 G_warning(_("Taking too long to find points for interpolation - " 00112 "please change the region to area where your points are. " 00113 "Continuing calculations...")); 00114 break; 00115 } 00116 i++; 00117 if (npt > params->KMAX2) 00118 /* decrease window */ 00119 { 00120 MAXENC = 1; 00121 nptprev = npt; 00122 temp1 = distxp; 00123 distxp = distx; 00124 distx = distxp - fabs(distx - temp1) * 0.5; 00125 temp2 = distyp; 00126 distyp = disty; 00127 disty = distyp - fabs(disty - temp2) * 0.5; 00128 /* decrease by 50% of a previous change in window */ 00129 } 00130 else { 00131 nptprev = npt; 00132 temp1 = distyp; 00133 distyp = disty; 00134 temp2 = distxp; 00135 distxp = distx; 00136 if (MAXENC) { 00137 disty = fabs(disty - temp1) * 0.5 + distyp; 00138 distx = fabs(distx - temp2) * 0.5 + distxp; 00139 } 00140 else { 00141 distx += distx; 00142 disty += disty; 00143 } 00144 /* decrease by 50% of extra distance */ 00145 } 00146 data->x_orig = xmn - distx; /* update window */ 00147 data->y_orig = ymn - disty; 00148 data->xmax = xmx + distx; 00149 data->ymax = ymx + disty; 00150 data->n_points = 0; 00151 npt = MT_region_data(info, info->root, data, params->KMAX2, 4); 00152 } 00153 00154 if (totsegm != 0) { 00155 G_percent(cursegm, totsegm, 1); 00156 } 00157 data->n_rows = ((struct quaddata *)(tree->data))->n_rows; 00158 data->n_cols = ((struct quaddata *)(tree->data))->n_cols; 00159 00160 /* for printing out overlapping segments */ 00161 ((struct quaddata *)(tree->data))->x_orig = xmn - distx; 00162 ((struct quaddata *)(tree->data))->y_orig = ymn - disty; 00163 ((struct quaddata *)(tree->data))->xmax = xmx + distx; 00164 ((struct quaddata *)(tree->data))->ymax = ymx + disty; 00165 00166 data->x_orig = xmn; 00167 data->y_orig = ymn; 00168 data->xmax = xmx; 00169 data->ymax = ymx; 00170 00171 if (!matrix) { 00172 if (! 00173 (matrix = 00174 G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) { 00175 G_warning(_("Out of memory")); 00176 return -1; 00177 } 00178 } 00179 if (!indx) { 00180 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) { 00181 G_warning(_("Out of memory")); 00182 return -1; 00183 } 00184 } 00185 if (!b) { 00186 if (!(b = G_alloc_vector(params->KMAX2 + 3))) { 00187 G_warning(_("Out of memory")); 00188 return -1; 00189 } 00190 } 00191 /* allocate memory for CV points only if cv is performed */ 00192 if (params->cv) { 00193 if (! 00194 (point = 00195 (struct triple *)G_malloc(sizeof(struct triple) * 00196 data->n_points))) { 00197 G_warning(_("Out of memory")); 00198 return -1; 00199 } 00200 } 00201 00202 /*normalize the data so that the side of average segment is about 1m */ 00203 /* put data_points into point only if CV is performed */ 00204 00205 for (i = 0; i < data->n_points; i++) { 00206 data->points[i].x = (data->points[i].x - data->x_orig) / dnorm; 00207 data->points[i].y = (data->points[i].y - data->y_orig) / dnorm; 00208 if (params->cv) { 00209 point[i].x = data->points[i].x; /*cv stuff */ 00210 point[i].y = data->points[i].y; /*cv stuff */ 00211 point[i].z = data->points[i].z; /*cv stuff */ 00212 } 00213 00214 /* commented out by Helena january 1997 as this is not necessary 00215 although it may be useful to put normalization of z back? 00216 data->points[i].z = data->points[i].z / dnorm; 00217 this made smoothing self-adjusting based on dnorm 00218 if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm; 00219 */ 00220 } 00221 00222 /* cv stuff */ 00223 if (params->cv) 00224 m_skip = data->n_points; 00225 else 00226 m_skip = 1; 00227 00228 /* remove after cleanup - this is just for testing */ 00229 skip_point.x = 0.; 00230 skip_point.y = 0.; 00231 skip_point.z = 0.; 00232 00233 for (skip_index = 0; skip_index < m_skip; skip_index++) { 00234 if (params->cv) { 00235 segtest = 0; 00236 j = 0; 00237 xx = point[skip_index].x * dnorm + data->x_orig + 00238 params->x_orig; 00239 yy = point[skip_index].y * dnorm + data->y_orig + 00240 params->y_orig; 00241 zz = point[skip_index].z; 00242 if (xx >= data->x_orig + params->x_orig && 00243 xx <= data->xmax + params->x_orig && 00244 yy >= data->y_orig + params->y_orig && 00245 yy <= data->ymax + params->y_orig) { 00246 segtest = 1; 00247 skip_point.x = point[skip_index].x; 00248 skip_point.y = point[skip_index].y; 00249 skip_point.z = point[skip_index].z; 00250 for (k = 0; k < m_skip; k++) { 00251 if (k != skip_index && params->cv) { 00252 data->points[j].x = point[k].x; 00253 data->points[j].y = point[k].y; 00254 data->points[j].z = point[k].z; 00255 j++; 00256 } 00257 } 00258 } /* segment area test */ 00259 } 00260 if (!params->cv) { 00261 if (params-> 00262 matrix_create(params, data->points, data->n_points, 00263 matrix, indx) < 0) 00264 return -1; 00265 } 00266 else if (segtest == 1) { 00267 if (params-> 00268 matrix_create(params, data->points, data->n_points - 1, 00269 matrix, indx) < 0) 00270 return -1; 00271 } 00272 if (!params->cv) { 00273 for (i = 0; i < data->n_points; i++) 00274 b[i + 1] = data->points[i].z; 00275 b[0] = 0.; 00276 G_lubksb(matrix, data->n_points + 1, indx, b); 00277 /* put here condition to skip ertot if not needed */ 00278 params->check_points(params, data, b, ertot, zmin, dnorm, 00279 skip_point); 00280 } 00281 else if (segtest == 1) { 00282 for (i = 0; i < data->n_points - 1; i++) 00283 b[i + 1] = data->points[i].z; 00284 b[0] = 0.; 00285 G_lubksb(matrix, data->n_points, indx, b); 00286 params->check_points(params, data, b, ertot, zmin, dnorm, 00287 skip_point); 00288 } 00289 } /*end of cv loop */ 00290 00291 if (!params->cv) 00292 if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) || 00293 (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) || 00294 (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) { 00295 00296 if (params->grid_calc(params, data, bitmask, 00297 zmin, zmax, zminac, zmaxac, gmin, gmax, 00298 c1min, c1max, c2min, c2max, ertot, b, 00299 offset1, dnorm) < 0) 00300 return -1; 00301 } 00302 00303 /* show after to catch 100% */ 00304 cursegm++; 00305 if (totsegm < cursegm) 00306 G_debug(1, "%d %d", totsegm, cursegm); 00307 00308 if (totsegm != 0) { 00309 G_percent(cursegm, totsegm, 1); 00310 } 00311 /* 00312 G_free_matrix(matrix); 00313 G_free_ivector(indx); 00314 G_free_vector(b); 00315 */ 00316 G_free(data->points); 00317 G_free(data); 00318 } 00319 return 1; 00320 } 00321 00322 static double smallest_segment(struct multtree *tree, int n_leafs) 00323 { 00324 static int first_time = 1; 00325 int ii; 00326 static double minside; 00327 double side; 00328 00329 if (tree == NULL) 00330 return 0; 00331 if (tree->data == NULL) 00332 return 0; 00333 if (tree->leafs != NULL) { 00334 for (ii = 0; ii < n_leafs; ii++) { 00335 side = smallest_segment(tree->leafs[ii], n_leafs); 00336 if (first_time) { 00337 minside = side; 00338 first_time = 0; 00339 } 00340 if (side < minside) 00341 minside = side; 00342 } 00343 } 00344 else { 00345 side = ((struct quaddata *)(tree->data))->xmax - 00346 ((struct quaddata *)(tree->data))->x_orig; 00347 return side; 00348 } 00349 00350 return minside; 00351 }