GRASS Programmer's Manual  6.4.2(2012)
dataquad.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines