GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * 00004 * Original program and various modifications: 00005 * Lubos Mitas 00006 * 00007 * GRASS4.1 version of the program and GRASS4.2 modifications: 00008 * H. Mitasova, 00009 * I. Kosinovsky, D. Gerdes 00010 * D. McCauley 00011 * 00012 * Copyright 1993, 1995: 00013 * L. Mitas , 00014 * H. Mitasova, 00015 * I. Kosinovsky, 00016 * D.Gerdes 00017 * D. McCauley 00018 * 00019 * modified by McCauley in August 1995 00020 * modified by Mitasova in August 1995, Nov. 1996 00021 * bug fixes(mask) and modif. for variable smoothing Mitasova Jan 1997 00022 * 00023 */ 00024 00025 00026 #include <stdio.h> 00027 #include <math.h> 00028 #include <unistd.h> 00029 #include <grass/gis.h> 00030 #include <grass/glocale.h> 00031 #include <grass/bitmap.h> 00032 00033 #include <grass/interpf.h> 00034 00035 00036 #define CEULER .57721566 00037 00038 00039 int IL_grid_calc_2d(struct interp_params *params, struct quaddata *data, /* given segment */ 00040 struct BM *bitmask, /* bitmask */ 00041 double zmin, double zmax, /* min and max input z-values */ 00042 double *zminac, double *zmaxac, /* min and max interp. z-values */ 00043 double *gmin, double *gmax, /* min and max inperp. slope val. */ 00044 double *c1min, double *c1max, double *c2min, double *c2max, /* min and max interp. curv. val. */ 00045 double *ertot, /* total interplating func. error */ 00046 double *b, /* solutions of linear equations */ 00047 int offset1, /* offset for temp file writing */ 00048 double dnorm) 00049 00050 /* 00051 * Calculates grid for the given segment represented by data (contains 00052 * n_rows, n_cols, ew_res,ns_res, and all points inside + overlap) using 00053 * solutions of system of lin. equations and interpolating functions 00054 * interp() and interpder(). Also calls secpar() to compute slope, aspect 00055 * and curvatures if required. 00056 */ 00057 { 00058 00059 /* 00060 * C C INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul. 00061 * c 00062 */ 00063 double x_or = data->x_orig; 00064 double y_or = data->y_orig; 00065 int n_rows = data->n_rows; 00066 int n_cols = data->n_cols; 00067 int n_points = data->n_points; 00068 struct triple *points; 00069 static double *w2 = NULL; 00070 static double *w = NULL; 00071 int cond1, cond2; 00072 double r; 00073 double stepix, stepiy, xx, xg, yg, xx2; 00074 double rfsta2, /* cons, cons1, */ wm, dx, dy, dxx, dyy, dxy, h, bmgd1, 00075 bmgd2; 00076 double r2, gd1, gd2; /* for interpder() */ 00077 int n1, k, l, m; 00078 int ngstc, nszc, ngstr, nszr; 00079 double zz; 00080 int bmask = 1; 00081 static int first_time_z = 1; 00082 int offset, offset2; 00083 double fstar2 = params->fi * params->fi / 4.; 00084 double tfsta2, tfstad; 00085 double ns_res, ew_res; 00086 double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */ 00087 double xxr, yyr; 00088 00089 if (params->theta) { 00090 teta = params->theta / 57.295779; /* deg to rad */ 00091 rsin = sin(teta); 00092 rcos = cos(teta); 00093 } 00094 if (params->scalex) 00095 scale = params->scalex; 00096 00097 ns_res = (((struct quaddata *)(data))->ymax - 00098 ((struct quaddata *)(data))->y_orig) / data->n_rows; 00099 ew_res = (((struct quaddata *)(data))->xmax - 00100 ((struct quaddata *)(data))->x_orig) / data->n_cols; 00101 00102 /* tfsta2 = fstar2 * 2.; modified after removing normalization of z */ 00103 tfsta2 = (fstar2 * 2.) / dnorm; 00104 tfstad = tfsta2 / dnorm; 00105 points = data->points; 00106 00107 /* 00108 * normalization 00109 */ 00110 stepix = ew_res / dnorm; 00111 stepiy = ns_res / dnorm; 00112 00113 cond2 = ((params->adxx != NULL) || (params->adyy != NULL) || 00114 (params->adxy != NULL)); 00115 cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2); 00116 00117 if (!w) { 00118 if (!(w = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) { 00119 G_warning(_("Out of memory")); 00120 return -1; 00121 } 00122 } 00123 if (!w2) { 00124 if (!(w2 = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) { 00125 G_warning(_("Out of memory")); 00126 return -1; 00127 } 00128 } 00129 n1 = n_points + 1; 00130 /* 00131 * C C INTERPOLATION * MOST INNER LOOPS ! C 00132 */ 00133 ngstc = (int)(x_or / ew_res + 0.5) + 1; 00134 nszc = ngstc + n_cols - 1; 00135 ngstr = (int)(y_or / ns_res + 0.5) + 1; 00136 nszr = ngstr + n_rows - 1; 00137 00138 00139 for (k = ngstr; k <= nszr; k++) { 00140 offset = offset1 * (k - 1); /* rows offset */ 00141 yg = (k - ngstr) * stepiy + stepiy / 2.; /* fixed by J.H. in July 01 */ 00142 for (m = 1; m <= n_points; m++) { 00143 wm = yg - points[m - 1].y; 00144 w[m] = wm; 00145 w2[m] = wm * wm; 00146 } 00147 for (l = ngstc; l <= nszc; l++) { 00148 if (bitmask != NULL) 00149 /* if(params->maskmap != NULL) PK Apr 03 MASK support */ 00150 bmask = BM_get(bitmask, l - 1, k - 1); /*fixed by helena jan 97 */ 00151 /* if(bmask==0 || bmask==-1) fprintf(stderr, "bmask=%d, at (%d,%d)\n", bmask, l, k); */ 00152 xg = (l - ngstc) * stepix + stepix / 2.; /*fixed by J.H. in July 01 */ 00153 dx = 0.; 00154 dy = 0.; 00155 dxx = 0.; 00156 dyy = 0.; 00157 dxy = 0.; 00158 zz = 0.; 00159 if (bmask == 1) { /* compute everything for area which is 00160 * not masked out */ 00161 h = b[0]; 00162 for (m = 1; m <= n_points; m++) { 00163 xx = xg - points[m - 1].x; 00164 if ((params->theta) && (params->scalex)) { 00165 /* we run anisotropy */ 00166 xxr = xx * rcos + w[m] * rsin; 00167 yyr = w[m] * rcos - xx * rsin; 00168 xx2 = xxr * xxr; 00169 w2[m] = yyr * yyr; 00170 r2 = scale * xx2 + w2[m]; 00171 r = r2; 00172 rfsta2 = scale * xx2 + w2[m]; 00173 } 00174 else { 00175 xx2 = xx * xx; 00176 r2 = xx2 + w2[m]; 00177 r = r2; 00178 rfsta2 = xx2 + w2[m]; 00179 } 00180 00181 h = h + b[m] * params->interp(r, params->fi); 00182 if (cond1) { 00183 if (!params->interpder(r, params->fi, &gd1, &gd2)) 00184 return -1; 00185 bmgd1 = b[m] * gd1; 00186 dx = dx + bmgd1 * xx; 00187 dy = dy + bmgd1 * w[m]; 00188 if (cond2) { 00189 bmgd2 = b[m] * gd2; 00190 dxx = dxx + bmgd2 * xx2 + bmgd1; 00191 dyy = dyy + bmgd2 * w2[m] + bmgd1; 00192 dxy = dxy + bmgd2 * xx * w[m]; 00193 } 00194 } 00195 } 00196 00197 /* zz = (h * dnorm) + zmin; replaced by helena jan. 97 due to removing norma 00198 lization of z and zm in segmen2d.c */ 00199 zz = h + zmin; 00200 if (first_time_z) { 00201 first_time_z = 0; 00202 *zmaxac = *zminac = zz; 00203 } 00204 *zmaxac = amax1(zz, *zmaxac); 00205 *zminac = amin1(zz, *zminac); 00206 if ((zz > zmax + 0.1 * (zmax - zmin)) 00207 || (zz < zmin - 0.1 * (zmax - zmin))) { 00208 static int once = 0; 00209 00210 if (!once) { 00211 once = 1; 00212 G_warning(_("Overshoot - increase in tension suggested. " 00213 "Overshoot occures at (%d,%d) cell. " 00214 "Z-value %f, zmin %f, zmax %f."), 00215 l, k, zz, zmin, zmax); 00216 } 00217 } 00218 00219 params->az[l] = (FCELL) zz; 00220 00221 if (cond1) { 00222 params->adx[l] = (FCELL) (-dx * tfsta2); 00223 params->ady[l] = (FCELL) (-dy * tfsta2); 00224 if (cond2) { 00225 params->adxx[l] = (FCELL) (-dxx * tfstad); 00226 params->adyy[l] = (FCELL) (-dyy * tfstad); 00227 params->adxy[l] = (FCELL) (-dxy * tfstad); 00228 } 00229 } 00230 00231 } 00232 else { 00233 G_set_d_null_value(params->az + l, 1); 00234 /* fprintf (stderr, "zz=%f, az[l]=%f, c=%d\n", zz, params->az[l], l); */ 00235 00236 if (cond1) { 00237 G_set_d_null_value(params->adx + l, 1); 00238 G_set_d_null_value(params->ady + l, 1); 00239 if (cond2) { 00240 G_set_d_null_value(params->adxx + l, 1); 00241 G_set_d_null_value(params->adyy + l, 1); 00242 G_set_d_null_value(params->adxy + l, 1); 00243 } 00244 } 00245 } 00246 00247 } 00248 if (cond1 && (params->deriv != 1)) { 00249 if (params->secpar(params, ngstc, nszc, k, bitmask, 00250 gmin, gmax, c1min, c1max, c2min, c2max, cond1, 00251 cond2) < 0) 00252 return -1; 00253 } 00254 00255 offset2 = (offset + ngstc - 1) * sizeof(FCELL); 00256 if (params->wr_temp(params, ngstc, nszc, offset2) < 0) 00257 return -1; 00258 } 00259 return 1; 00260 }