GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 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 H.Mitasova November 1996 to include variable smoothing 00010 * 00011 */ 00012 00013 #include <stdio.h> 00014 #include <stdlib.h> 00015 #include <grass/dataquad.h> 00016 00017 /* sm added to point structure */ 00018 struct triple *quad_point_new(double x, double y, double z, double sm) 00019 /* Initializes POINT structure with given arguments */ 00020 { 00021 struct triple *point; 00022 00023 if (!(point = (struct triple *)malloc(sizeof(struct triple)))) { 00024 return NULL; 00025 } 00026 00027 point->x = x; 00028 point->y = y; 00029 point->z = z; 00030 point->sm = sm; 00031 00032 return point; 00033 } 00034 00035 00036 struct quaddata *quad_data_new(double x_or, double y_or, double xmax, 00037 double ymax, int rows, int cols, int n_points, 00038 int kmax) 00039 /* Initializes QUADDATA structure with given arguments */ 00040 { 00041 struct quaddata *data; 00042 int i; 00043 00044 if (!(data = (struct quaddata *)malloc(sizeof(struct quaddata)))) { 00045 return NULL; 00046 } 00047 00048 data->x_orig = x_or; 00049 data->y_orig = y_or; 00050 data->xmax = xmax; 00051 data->ymax = ymax; 00052 data->n_rows = rows; 00053 data->n_cols = cols; 00054 data->n_points = n_points; 00055 data->points = 00056 (struct triple *)malloc(sizeof(struct triple) * (kmax + 1)); 00057 if (!data->points) 00058 return NULL; 00059 for (i = 0; i <= kmax; i++) { 00060 data->points[i].x = 0.; 00061 data->points[i].y = 0.; 00062 data->points[i].z = 0.; 00063 data->points[i].sm = 0.; 00064 } 00065 00066 return data; 00067 } 00068 00069 00070 00071 00072 00073 int quad_compare(struct triple *point, struct quaddata *data) 00074 /* returns the quadrant the point should be inserted in */ 00075 /* called by divide() */ 00076 { 00077 int cond1, cond2, cond3, cond4, rows, cols; 00078 double ew_res, ns_res; 00079 00080 ew_res = (data->xmax - data->x_orig) / data->n_cols; 00081 ns_res = (data->ymax - data->y_orig) / data->n_rows; 00082 00083 00084 if (data == NULL) 00085 return -1; 00086 if (data->n_rows % 2 == 0) { 00087 rows = data->n_rows / 2; 00088 } 00089 else { 00090 rows = (int)(data->n_rows / 2) + 1; 00091 } 00092 00093 if (data->n_cols % 2 == 0) { 00094 cols = data->n_cols / 2; 00095 } 00096 else { 00097 cols = (int)(data->n_cols / 2) + 1; 00098 } 00099 cond1 = (point->x >= data->x_orig); 00100 cond2 = (point->x >= data->x_orig + ew_res * cols); 00101 cond3 = (point->y >= data->y_orig); 00102 cond4 = (point->y >= data->y_orig + ns_res * rows); 00103 if (cond1 && cond3) { 00104 if (cond2 && cond4) 00105 return NE; 00106 if (cond2) 00107 return SE; 00108 if (cond4) 00109 return NW; 00110 return SW; 00111 } 00112 else 00113 return 0; 00114 } 00115 00116 00117 int quad_add_data(struct triple *point, struct quaddata *data, double dmin) 00118 /* Adds POINT to a given DATA . Called by tree function insert_quad() */ 00119 /* and by data function quad_divide_data() */ 00120 { 00121 int n, i, cond; 00122 double xx, yy, r; 00123 00124 cond = 1; 00125 if (data == NULL) { 00126 fprintf(stderr, "add_data: data is NULL \n"); 00127 return -5; 00128 } 00129 for (i = 0; i < data->n_points; i++) { 00130 xx = data->points[i].x - point->x; 00131 yy = data->points[i].y - point->y; 00132 r = xx * xx + yy * yy; 00133 if (r <= dmin) { 00134 cond = 0; 00135 break; 00136 } 00137 } 00138 00139 if (cond) { 00140 n = (data->n_points)++; 00141 data->points[n].x = point->x; 00142 data->points[n].y = point->y; 00143 data->points[n].z = point->z; 00144 data->points[n].sm = point->sm; 00145 } 00146 return cond; 00147 } 00148 00149 00150 00151 00152 int quad_intersect(struct quaddata *data_inter, struct quaddata *data) 00153 /* Checks if region defined by DATA intersects the region defined 00154 by data_inter. Called by tree function MT_region_data() */ 00155 { 00156 double xmin, xmax, ymin, ymax; 00157 00158 xmin = data_inter->x_orig; 00159 xmax = data_inter->xmax; 00160 ymin = data_inter->y_orig; 00161 ymax = data_inter->ymax; 00162 00163 if (((data->x_orig >= xmin) && (data->x_orig <= xmax) 00164 && (((data->y_orig >= ymin) && (data->y_orig <= ymax)) 00165 || ((ymin >= data->y_orig) && (ymin <= data->ymax)) 00166 ) 00167 ) 00168 || ((xmin >= data->x_orig) && (xmin <= data->xmax) 00169 && (((ymin >= data->y_orig) && (ymin <= data->ymax)) 00170 || ((data->y_orig >= ymin) && (data->y_orig <= ymax)) 00171 ) 00172 ) 00173 ) { 00174 return 1; 00175 } 00176 else 00177 return 0; 00178 } 00179 00180 00181 00182 00183 int quad_division_check(struct quaddata *data, int kmax) 00184 /* Checks if DATA needs to be divided. If data->points is empty, 00185 returns -1; if its not empty but there aren't enough points 00186 in DATA for division returns 0. Othervise (if its not empty and 00187 there are too many points) returns 1. Called by MT_insert() */ 00188 { 00189 if (data->points == NULL) 00190 return -1; 00191 if (data->n_points < kmax) 00192 return 0; 00193 else 00194 return 1; 00195 } 00196 00197 00198 00199 struct quaddata **quad_divide_data(struct quaddata *data, int kmax, 00200 double dmin) 00201 /* Divides DATA into 4 new datas reinserting data->points in 00202 them by calling data function quad_compare() to detrmine 00203 were to insert. Called by MT_divide(). Returns array of 4 new datas */ 00204 { 00205 struct quaddata **datas; 00206 int cols1, cols2, rows1, rows2, i; /*j1, j2, jmin = 0; */ 00207 double dx, dy; /* x2, y2, dist, mindist; */ 00208 double xr, xm, xl, yr, ym, yl; /* left, right, middle coord */ 00209 double ew_res, ns_res; 00210 00211 ew_res = (data->xmax - data->x_orig) / data->n_cols; 00212 ns_res = (data->ymax - data->y_orig) / data->n_rows; 00213 00214 if ((data->n_cols <= 1) || (data->n_rows <= 1)) { 00215 fprintf(stderr, 00216 "Points are too concentrated -- please increase DMIN\n"); 00217 exit(0); 00218 } 00219 00220 if (data->n_cols % 2 == 0) { 00221 cols1 = data->n_cols / 2; 00222 cols2 = cols1; 00223 } 00224 else { 00225 cols2 = (int)(data->n_cols / 2); 00226 cols1 = cols2 + 1; 00227 } 00228 if (data->n_rows % 2 == 0) { 00229 rows1 = data->n_rows / 2; 00230 rows2 = rows1; 00231 } 00232 else { 00233 rows2 = (int)(data->n_rows / 2); 00234 rows1 = rows2 + 1; 00235 } 00236 00237 dx = cols1 * ew_res; 00238 dy = rows1 * ns_res; 00239 00240 xl = data->x_orig; 00241 xm = xl + dx; 00242 xr = data->xmax; 00243 yl = data->y_orig; 00244 ym = yl + dy; 00245 yr = data->ymax; 00246 00247 if (!(datas = (struct quaddata **)malloc(sizeof(struct quaddata *) * 5))) { 00248 return NULL; 00249 } 00250 datas[NE] = quad_data_new(xm, ym, xr, yr, rows2, cols2, 0, kmax); 00251 datas[SW] = quad_data_new(xl, yl, xm, ym, rows1, cols1, 0, kmax); 00252 datas[SE] = quad_data_new(xm, yl, xr, ym, rows1, cols2, 0, kmax); 00253 datas[NW] = quad_data_new(xl, ym, xm, yr, rows2, cols1, 0, kmax); 00254 for (i = 0; i < data->n_points; i++) { 00255 switch (quad_compare(data->points + i, data)) { 00256 case SW: 00257 { 00258 quad_add_data(data->points + i, datas[SW], dmin); 00259 break; 00260 } 00261 case SE: 00262 { 00263 quad_add_data(data->points + i, datas[SE], dmin); 00264 break; 00265 } 00266 case NW: 00267 { 00268 quad_add_data(data->points + i, datas[NW], dmin); 00269 break; 00270 } 00271 case NE: 00272 { 00273 quad_add_data(data->points + i, datas[NE], dmin); 00274 break; 00275 } 00276 } 00277 } 00278 data->points = NULL; 00279 return datas; 00280 } 00281 00282 00283 00284 00285 int 00286 quad_get_points(struct quaddata *data_inter, struct quaddata *data, int MAX) 00287 /* Gets such points from DATA that lie within region determined by 00288 data_inter. Called by tree function region_data(). */ 00289 { 00290 int i, ind; 00291 int n = 0; 00292 int l = 0; 00293 double xmin, xmax, ymin, ymax; 00294 struct triple *point; 00295 00296 xmin = data_inter->x_orig; 00297 xmax = data_inter->xmax; 00298 ymin = data_inter->y_orig; 00299 ymax = data_inter->ymax; 00300 for (i = 0; i < data->n_points; i++) { 00301 point = data->points + i; 00302 if (l >= MAX) 00303 return MAX + 1; 00304 if ((point->x > xmin) && (point->x < xmax) 00305 && (point->y > ymin) && (point->y < ymax)) { 00306 ind = data_inter->n_points++; 00307 data_inter->points[ind].x = point->x; 00308 data_inter->points[ind].y = point->y; 00309 data_inter->points[ind].z = point->z; 00310 data_inter->points[ind].sm = point->sm; 00311 l = l + 1; 00312 00313 } 00314 } 00315 n = l; 00316 return (n); 00317 }