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 * 00022 */ 00023 00024 00025 #include <stdio.h> 00026 #include <math.h> 00027 #include <unistd.h> 00028 #include <grass/gis.h> 00029 #include <grass/site.h> 00030 #include <grass/Vect.h> 00031 #include <grass/dbmi.h> 00032 00033 #define POINT2D_C 00034 #include <grass/interpf.h> 00035 00036 int IL_check_at_points_2d(struct interp_params *params, struct quaddata *data, /* current region */ 00037 double *b, /* solution of linear equations */ 00038 double *ertot, /* total error */ 00039 double zmin, /* min z-value */ 00040 double dnorm, struct triple skip_point) 00041 00042 /* 00043 * Checks if interpolating function interp() evaluates correct z-values at 00044 * given points. If smoothing is used calculate the maximum error caused 00045 * by smoothing. 00046 */ 00047 { 00048 int n_points = data->n_points; /* number of points */ 00049 struct triple *points = data->points; /* points for interpolation */ 00050 double east = data->xmax; 00051 double west = data->x_orig; 00052 double north = data->ymax; 00053 double south = data->y_orig; 00054 double rfsta2, errmax, h, xx, yy, r2, hz, zz, err, xmm, ymm, r; 00055 double skip_err; 00056 int n1, mm, m, cat; 00057 double fstar2; 00058 int inside; 00059 00060 /* Site *site; */ 00061 char buf[1024]; 00062 00063 00064 /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL) 00065 G_fatal_error ("Memory error for site struct"); */ 00066 00067 fstar2 = params->fi * params->fi / 4.; 00068 errmax = .0; 00069 n1 = n_points + 1; 00070 for (mm = 1; mm <= n_points; mm++) { 00071 h = b[0]; 00072 for (m = 1; m <= n_points; m++) { 00073 xx = points[mm - 1].x - points[m - 1].x; 00074 yy = points[mm - 1].y - points[m - 1].y; 00075 r2 = yy * yy + xx * xx; 00076 if (r2 != 0.) { 00077 rfsta2 = fstar2 * r2; 00078 r = r2; 00079 h = h + b[m] * params->interp(r, params->fi); 00080 } 00081 } 00082 /* modified by helena january 1997 - normalization of z was 00083 removed from segm2d.c and interp2d.c 00084 hz = (h * dnorm) + zmin; 00085 zz = (points[mm - 1].z * dnorm) + zmin; 00086 */ 00087 hz = h + zmin; 00088 zz = points[mm - 1].z + zmin; 00089 err = hz - zz; 00090 xmm = points[mm - 1].x * dnorm + params->x_orig + west; 00091 ymm = points[mm - 1].y * dnorm + params->y_orig + south; 00092 if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) 00093 && (ymm >= south + params->y_orig) && 00094 (ymm <= north + params->y_orig)) 00095 inside = 1; 00096 else 00097 inside = 0; 00098 00099 if (params->fddevi != NULL) { 00100 00101 if (inside) { /* if the point is inside the region */ 00102 Vect_reset_line(Pnts); 00103 Vect_reset_cats(Cats2); 00104 00105 Vect_append_point(Pnts, xmm, ymm, zz); 00106 cat = count; 00107 Vect_cat_set(Cats2, 1, cat); 00108 Vect_write_line(&Map2, GV_POINT, Pnts, Cats2); 00109 00110 db_zero_string(&sql2); 00111 sprintf(buf, "insert into %s values ( %d ", ff->table, cat); 00112 db_append_string(&sql2, buf); 00113 00114 sprintf(buf, ", %f", err); 00115 db_append_string(&sql2, buf); 00116 db_append_string(&sql2, ")"); 00117 G_debug(3, db_get_string(&sql2)); 00118 00119 if (db_execute_immediate(driver2, &sql2) != DB_OK) { 00120 db_close_database(driver2); 00121 db_shutdown_driver(driver2); 00122 G_fatal_error("Cannot insert new row: %s", 00123 db_get_string(&sql2)); 00124 } 00125 count++; 00126 00127 } 00128 } 00129 (*ertot) += err * err; 00130 } 00131 00132 /* cv stuff */ 00133 if (params->cv) { 00134 h = b[0]; 00135 for (m = 1; m <= n_points - 1; m++) { 00136 xx = points[m - 1].x - skip_point.x; 00137 yy = points[m - 1].y - skip_point.y; 00138 r2 = yy * yy + xx * xx; 00139 if (r2 != 0.) { 00140 rfsta2 = fstar2 * r2; 00141 r = r2; 00142 h = h + b[m] * params->interp(r, params->fi); 00143 } 00144 } 00145 hz = h + zmin; 00146 zz = skip_point.z + zmin; 00147 skip_err = hz - zz; 00148 xmm = skip_point.x * dnorm + params->x_orig + west; 00149 ymm = skip_point.y * dnorm + params->y_orig + south; 00150 00151 if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) 00152 && (ymm >= south + params->y_orig) && 00153 (ymm <= north + params->y_orig)) 00154 inside = 1; 00155 else 00156 inside = 0; 00157 00158 if (inside) { /* if the point is inside the region */ 00159 Vect_reset_line(Pnts); 00160 Vect_reset_cats(Cats2); 00161 00162 Vect_append_point(Pnts, xmm, ymm, zz); 00163 cat = count; 00164 Vect_cat_set(Cats2, 1, cat); 00165 Vect_write_line(&Map2, GV_POINT, Pnts, Cats2); 00166 00167 db_zero_string(&sql2); 00168 sprintf(buf, "insert into %s values ( %d ", ff->table, cat); 00169 db_append_string(&sql2, buf); 00170 00171 sprintf(buf, ", %f", skip_err); 00172 db_append_string(&sql2, buf); 00173 db_append_string(&sql2, ")"); 00174 G_debug(3, db_get_string(&sql2)); 00175 00176 if (db_execute_immediate(driver2, &sql2) != DB_OK) { 00177 db_close_database(driver2); 00178 db_shutdown_driver(driver2); 00179 G_fatal_error("Cannot insert new row: %s", 00180 db_get_string(&sql2)); 00181 } 00182 count++; 00183 } 00184 } /* cv */ 00185 00186 00187 return 1; 00188 }