GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
dataquad.c
Go to the documentation of this file.
1 
2 /*-
3  * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 1993
4  * University of Illinois
5  * US Army Construction Engineering Research Lab
6  * Copyright 1993, H. Mitasova (University of Illinois),
7  * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
8  *
9  * Modified by H.Mitasova November 1996 to include variable smoothing
10  *
11  */
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <grass/dataquad.h>
16 
17 /* sm added to point structure */
18 struct triple *quad_point_new(double x, double y, double z, double sm)
19 /* Initializes POINT structure with given arguments */
20 {
21  struct triple *point;
22 
23  if (!(point = (struct triple *)malloc(sizeof(struct triple)))) {
24  return NULL;
25  }
26 
27  point->x = x;
28  point->y = y;
29  point->z = z;
30  point->sm = sm;
31 
32  return point;
33 }
34 
35 
36 struct quaddata *quad_data_new(double x_or, double y_or, double xmax,
37  double ymax, int rows, int cols, int n_points,
38  int kmax)
39 /* Initializes QUADDATA structure with given arguments */
40 {
41  struct quaddata *data;
42  int i;
43 
44  if (!(data = (struct quaddata *)malloc(sizeof(struct quaddata)))) {
45  return NULL;
46  }
47 
48  data->x_orig = x_or;
49  data->y_orig = y_or;
50  data->xmax = xmax;
51  data->ymax = ymax;
52  data->n_rows = rows;
53  data->n_cols = cols;
54  data->n_points = n_points;
55  data->points =
56  (struct triple *)malloc(sizeof(struct triple) * (kmax + 1));
57  if (!data->points)
58  return NULL;
59  for (i = 0; i <= kmax; i++) {
60  data->points[i].x = 0.;
61  data->points[i].y = 0.;
62  data->points[i].z = 0.;
63  data->points[i].sm = 0.;
64  }
65 
66  return data;
67 }
68 
69 
70 
71 
72 
73 int quad_compare(struct triple *point, struct quaddata *data)
74 /* returns the quadrant the point should be inserted in */
75 /* called by divide() */
76 {
77  int cond1, cond2, cond3, cond4, rows, cols;
78  double ew_res, ns_res;
79 
80  ew_res = (data->xmax - data->x_orig) / data->n_cols;
81  ns_res = (data->ymax - data->y_orig) / data->n_rows;
82 
83 
84  if (data == NULL)
85  return -1;
86  if (data->n_rows % 2 == 0) {
87  rows = data->n_rows / 2;
88  }
89  else {
90  rows = (int)(data->n_rows / 2) + 1;
91  }
92 
93  if (data->n_cols % 2 == 0) {
94  cols = data->n_cols / 2;
95  }
96  else {
97  cols = (int)(data->n_cols / 2) + 1;
98  }
99  cond1 = (point->x >= data->x_orig);
100  cond2 = (point->x >= data->x_orig + ew_res * cols);
101  cond3 = (point->y >= data->y_orig);
102  cond4 = (point->y >= data->y_orig + ns_res * rows);
103  if (cond1 && cond3) {
104  if (cond2 && cond4)
105  return NE;
106  if (cond2)
107  return SE;
108  if (cond4)
109  return NW;
110  return SW;
111  }
112  else
113  return 0;
114 }
115 
116 
117 int quad_add_data(struct triple *point, struct quaddata *data, double dmin)
118 /* Adds POINT to a given DATA . Called by tree function insert_quad() */
119 /* and by data function quad_divide_data() */
120 {
121  int n, i, cond;
122  double xx, yy, r;
123 
124  cond = 1;
125  if (data == NULL) {
126  fprintf(stderr, "add_data: data is NULL \n");
127  return -5;
128  }
129  for (i = 0; i < data->n_points; i++) {
130  xx = data->points[i].x - point->x;
131  yy = data->points[i].y - point->y;
132  r = xx * xx + yy * yy;
133  if (r <= dmin) {
134  cond = 0;
135  break;
136  }
137  }
138 
139  if (cond) {
140  n = (data->n_points)++;
141  data->points[n].x = point->x;
142  data->points[n].y = point->y;
143  data->points[n].z = point->z;
144  data->points[n].sm = point->sm;
145  }
146  return cond;
147 }
148 
149 
150 
151 
152 int quad_intersect(struct quaddata *data_inter, struct quaddata *data)
153 /* Checks if region defined by DATA intersects the region defined
154  by data_inter. Called by tree function MT_region_data() */
155 {
156  double xmin, xmax, ymin, ymax;
157 
158  xmin = data_inter->x_orig;
159  xmax = data_inter->xmax;
160  ymin = data_inter->y_orig;
161  ymax = data_inter->ymax;
162 
163  if (((data->x_orig >= xmin) && (data->x_orig <= xmax)
164  && (((data->y_orig >= ymin) && (data->y_orig <= ymax))
165  || ((ymin >= data->y_orig) && (ymin <= data->ymax))
166  )
167  )
168  || ((xmin >= data->x_orig) && (xmin <= data->xmax)
169  && (((ymin >= data->y_orig) && (ymin <= data->ymax))
170  || ((data->y_orig >= ymin) && (data->y_orig <= ymax))
171  )
172  )
173  ) {
174  return 1;
175  }
176  else
177  return 0;
178 }
179 
180 
181 
182 
183 int quad_division_check(struct quaddata *data, int kmax)
184 /* Checks if DATA needs to be divided. If data->points is empty,
185  returns -1; if its not empty but there aren't enough points
186  in DATA for division returns 0. Othervise (if its not empty and
187  there are too many points) returns 1. Called by MT_insert() */
188 {
189  if (data->points == NULL)
190  return -1;
191  if (data->n_points < kmax)
192  return 0;
193  else
194  return 1;
195 }
196 
197 
198 
199 struct quaddata **quad_divide_data(struct quaddata *data, int kmax,
200  double dmin)
201 /* Divides DATA into 4 new datas reinserting data->points in
202  them by calling data function quad_compare() to detrmine
203  were to insert. Called by MT_divide(). Returns array of 4 new datas */
204 {
205  struct quaddata **datas;
206  int cols1, cols2, rows1, rows2, i; /*j1, j2, jmin = 0; */
207  double dx, dy; /* x2, y2, dist, mindist; */
208  double xr, xm, xl, yr, ym, yl; /* left, right, middle coord */
209  double ew_res, ns_res;
210 
211  ew_res = (data->xmax - data->x_orig) / data->n_cols;
212  ns_res = (data->ymax - data->y_orig) / data->n_rows;
213 
214  if ((data->n_cols <= 1) || (data->n_rows <= 1)) {
215  fprintf(stderr,
216  "Points are too concentrated -- please increase DMIN\n");
217  exit(0);
218  }
219 
220  if (data->n_cols % 2 == 0) {
221  cols1 = data->n_cols / 2;
222  cols2 = cols1;
223  }
224  else {
225  cols2 = (int)(data->n_cols / 2);
226  cols1 = cols2 + 1;
227  }
228  if (data->n_rows % 2 == 0) {
229  rows1 = data->n_rows / 2;
230  rows2 = rows1;
231  }
232  else {
233  rows2 = (int)(data->n_rows / 2);
234  rows1 = rows2 + 1;
235  }
236 
237  dx = cols1 * ew_res;
238  dy = rows1 * ns_res;
239 
240  xl = data->x_orig;
241  xm = xl + dx;
242  xr = data->xmax;
243  yl = data->y_orig;
244  ym = yl + dy;
245  yr = data->ymax;
246 
247  if (!(datas = (struct quaddata **)malloc(sizeof(struct quaddata *) * 5))) {
248  return NULL;
249  }
250  datas[NE] = quad_data_new(xm, ym, xr, yr, rows2, cols2, 0, kmax);
251  datas[SW] = quad_data_new(xl, yl, xm, ym, rows1, cols1, 0, kmax);
252  datas[SE] = quad_data_new(xm, yl, xr, ym, rows1, cols2, 0, kmax);
253  datas[NW] = quad_data_new(xl, ym, xm, yr, rows2, cols1, 0, kmax);
254  for (i = 0; i < data->n_points; i++) {
255  switch (quad_compare(data->points + i, data)) {
256  case SW:
257  {
258  quad_add_data(data->points + i, datas[SW], dmin);
259  break;
260  }
261  case SE:
262  {
263  quad_add_data(data->points + i, datas[SE], dmin);
264  break;
265  }
266  case NW:
267  {
268  quad_add_data(data->points + i, datas[NW], dmin);
269  break;
270  }
271  case NE:
272  {
273  quad_add_data(data->points + i, datas[NE], dmin);
274  break;
275  }
276  }
277  }
278  data->points = NULL;
279  return datas;
280 }
281 
282 
283 
284 
285 int
286 quad_get_points(struct quaddata *data_inter, struct quaddata *data, int MAX)
287 /* Gets such points from DATA that lie within region determined by
288  data_inter. Called by tree function region_data(). */
289 {
290  int i, ind;
291  int n = 0;
292  int l = 0;
293  double xmin, xmax, ymin, ymax;
294  struct triple *point;
295 
296  xmin = data_inter->x_orig;
297  xmax = data_inter->xmax;
298  ymin = data_inter->y_orig;
299  ymax = data_inter->ymax;
300  for (i = 0; i < data->n_points; i++) {
301  point = data->points + i;
302  if (l >= MAX)
303  return MAX + 1;
304  if ((point->x > xmin) && (point->x < xmax)
305  && (point->y > ymin) && (point->y < ymax)) {
306  ind = data_inter->n_points++;
307  data_inter->points[ind].x = point->x;
308  data_inter->points[ind].y = point->y;
309  data_inter->points[ind].z = point->z;
310  data_inter->points[ind].sm = point->sm;
311  l = l + 1;
312 
313  }
314  }
315  n = l;
316  return (n);
317 }
int l
Definition: dataquad.c:292
double xmax
Definition: dataquad.c:293
double y_orig
Definition: dataquad.h:33
struct triple * point
Definition: dataquad.c:294
int quad_intersect(struct quaddata *data_inter, struct quaddata *data)
Definition: dataquad.c:152
float r
Definition: named_colr.c:8
double xmin
Definition: dataquad.c:293
double z
Definition: dataquad.h:26
#define NE
Definition: dataquad.h:18
double ymin
Definition: dataquad.c:293
#define NW
Definition: dataquad.h:17
int y
Definition: plot.c:34
double x_orig
Definition: dataquad.h:32
DCELL dmin
Definition: g3dcolor.c:53
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition: dataquad.c:18
struct triple * points
Definition: dataquad.h:39
#define SW
Definition: dataquad.h:19
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
Definition: dataquad.c:36
int quad_compare(struct triple *point, struct quaddata *data)
Definition: dataquad.c:73
tuple data
void * malloc(YYSIZE_T)
struct quaddata ** quad_divide_data(struct quaddata *data, int kmax, double dmin)
Definition: dataquad.c:199
double x
Definition: dataquad.h:24
int
Definition: g3dcolor.c:48
double sm
Definition: dataquad.h:27
double ymax
Definition: dataquad.c:293
double ymax
Definition: dataquad.h:35
int quad_division_check(struct quaddata *data, int kmax)
Definition: dataquad.c:183
#define SE
Definition: dataquad.h:20
return NULL
Definition: dbfopen.c:1394
tuple cols
#define MAX(a, b)
double xmax
Definition: dataquad.h:34
int n_rows
Definition: dataquad.h:36
int n_cols
Definition: dataquad.h:37
int n
Definition: dataquad.c:291
int quad_add_data(struct triple *point, struct quaddata *data, double dmin)
Definition: dataquad.c:117
int quad_get_points(struct quaddata *, struct quaddata *, int)
double y
Definition: dataquad.h:25
int n_points
Definition: dataquad.h:38