GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
dataquad.c
Go to the documentation of this file.
1 /*!
2  * \file qtree.c
3  *
4  * \author
5  * H. Mitasova, I. Kosinovsky, D. Gerdes, Fall 1993,
6  * University of Illinois and
7  * US Army Construction Engineering Research Lab
8  *
9  * \author H. Mitasova (University of Illinois),
10  * \author I. Kosinovsky, (USA-CERL)
11  * \author D.Gerdes (USA-CERL)
12  *
13  * \author modified by H. Mitasova, November 1996 (include variable smoothing)
14  *
15  * \copyright
16  * (C) 1993-1996 by Helena Mitasova and the GRASS Development Team
17  *
18  * \copyright
19  * This program is free software under the
20  * GNU General Public License (>=v2).
21  * Read the file COPYING that comes with GRASS for details.
22  */
23 
24 
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <grass/dataquad.h>
28 
29 
30 /*!
31  * Initialize point structure with given arguments
32  *
33  * This is a constructor of the point structure and it allocates memory.
34  *
35  * \note
36  * Smoothing is part of the point structure
37  */
38 struct triple *quad_point_new(double x, double y, double z, double sm)
39 {
40  struct triple *point;
41 
42  if (!(point = (struct triple *)malloc(sizeof(struct triple)))) {
43  return NULL;
44  }
45 
46  point->x = x;
47  point->y = y;
48  point->z = z;
49  point->sm = sm;
50 
51  return point;
52 }
53 
54 
55 /*!
56  * Initialize quaddata structure with given arguments
57  *
58  * This is a constructor of the quaddata structure and it allocates memory.
59  * It also creates (and allocates memory for) the given number of points
60  * (given by *kmax*). The point attributes are set to zero.
61  */
62 struct quaddata *quad_data_new(double x_or, double y_or, double xmax,
63  double ymax, int rows, int cols, int n_points,
64  int kmax)
65 {
66  struct quaddata *data;
67  int i;
68 
69  if (!(data = (struct quaddata *)malloc(sizeof(struct quaddata)))) {
70  return NULL;
71  }
72 
73  data->x_orig = x_or;
74  data->y_orig = y_or;
75  data->xmax = xmax;
76  data->ymax = ymax;
77  data->n_rows = rows;
78  data->n_cols = cols;
79  data->n_points = n_points;
80  data->points =
81  (struct triple *)malloc(sizeof(struct triple) * (kmax + 1));
82  if (!data->points)
83  return NULL;
84  for (i = 0; i <= kmax; i++) {
85  data->points[i].x = 0.;
86  data->points[i].y = 0.;
87  data->points[i].z = 0.;
88  data->points[i].sm = 0.;
89  }
90 
91  return data;
92 }
93 
94 
95 /*!
96  * Return the quadrant the point should be inserted in
97  */
98 int quad_compare(struct triple *point, struct quaddata *data)
99 {
100  int cond1, cond2, cond3, cond4, rows, cols;
101  double ew_res, ns_res;
102 
103  ew_res = (data->xmax - data->x_orig) / data->n_cols;
104  ns_res = (data->ymax - data->y_orig) / data->n_rows;
105 
106 
107  if (data == NULL)
108  return -1;
109  if (data->n_rows % 2 == 0) {
110  rows = data->n_rows / 2;
111  }
112  else {
113  rows = (int)(data->n_rows / 2) + 1;
114  }
115 
116  if (data->n_cols % 2 == 0) {
117  cols = data->n_cols / 2;
118  }
119  else {
120  cols = (int)(data->n_cols / 2) + 1;
121  }
122  cond1 = (point->x >= data->x_orig);
123  cond2 = (point->x >= data->x_orig + ew_res * cols);
124  cond3 = (point->y >= data->y_orig);
125  cond4 = (point->y >= data->y_orig + ns_res * rows);
126  if (cond1 && cond3) {
127  if (cond2 && cond4)
128  return NE;
129  if (cond2)
130  return SE;
131  if (cond4)
132  return NW;
133  return SW;
134  }
135  else
136  return 0;
137 }
138 
139 
140 /*!
141  * Add point to a given *data*.
142  */
143 int quad_add_data(struct triple *point, struct quaddata *data, double dmin)
144 {
145  int n, i, cond;
146  double xx, yy, r;
147 
148  cond = 1;
149  if (data == NULL) {
150  fprintf(stderr, "add_data: data is NULL \n");
151  return -5;
152  }
153  for (i = 0; i < data->n_points; i++) {
154  xx = data->points[i].x - point->x;
155  yy = data->points[i].y - point->y;
156  r = xx * xx + yy * yy;
157  if (r <= dmin) {
158  cond = 0;
159  break;
160  }
161  }
162 
163  if (cond) {
164  n = (data->n_points)++;
165  data->points[n].x = point->x;
166  data->points[n].y = point->y;
167  data->points[n].z = point->z;
168  data->points[n].sm = point->sm;
169  }
170  return cond;
171 }
172 
173 
174 /*!
175  * Check intersection of two quaddata structures
176  *
177  * Checks if region defined by *data* intersects the region defined
178  * by *data_inter*.
179  */
180 int quad_intersect(struct quaddata *data_inter, struct quaddata *data)
181 {
182  double xmin, xmax, ymin, ymax;
183 
184  xmin = data_inter->x_orig;
185  xmax = data_inter->xmax;
186  ymin = data_inter->y_orig;
187  ymax = data_inter->ymax;
188 
189  if (((data->x_orig >= xmin) && (data->x_orig <= xmax)
190  && (((data->y_orig >= ymin) && (data->y_orig <= ymax))
191  || ((ymin >= data->y_orig) && (ymin <= data->ymax))
192  )
193  )
194  || ((xmin >= data->x_orig) && (xmin <= data->xmax)
195  && (((ymin >= data->y_orig) && (ymin <= data->ymax))
196  || ((data->y_orig >= ymin) && (data->y_orig <= ymax))
197  )
198  )
199  ) {
200  return 1;
201  }
202  else
203  return 0;
204 }
205 
206 
207 /*!
208  * Check if *data* needs to be divided
209  *
210  * Checks if *data* needs to be divided. If `data->points` is empty,
211  * returns -1; if its not empty but there aren't enough points
212  * in *data* for division returns 0. Otherwise (if its not empty and
213  * there are too many points) returns 1.
214  *
215  * \returns 1 if division is needed
216  * \returns 0 if division is not needed
217  * \returns -1 if there are no points
218  */
219 int quad_division_check(struct quaddata *data, int kmax)
220 {
221  if (data->points == NULL)
222  return -1;
223  if (data->n_points < kmax)
224  return 0;
225  else
226  return 1;
227 }
228 
229 
230 /*!
231  * Divide *data* into four new ones
232  *
233  * Divides *data* into 4 new datas reinserting `data->points` in
234  * them by calling data function `quad_compare()` to determine
235  * were to insert. Returns array of 4 new datas (allocates memory).
236  */
237 struct quaddata **quad_divide_data(struct quaddata *data, int kmax,
238  double dmin)
239 {
240  struct quaddata **datas;
241  int cols1, cols2, rows1, rows2, i; /*j1, j2, jmin = 0; */
242  double dx, dy; /* x2, y2, dist, mindist; */
243  double xr, xm, xl, yr, ym, yl; /* left, right, middle coord */
244  double ew_res, ns_res;
245 
246  ew_res = (data->xmax - data->x_orig) / data->n_cols;
247  ns_res = (data->ymax - data->y_orig) / data->n_rows;
248 
249  if ((data->n_cols <= 1) || (data->n_rows <= 1)) {
250  fprintf(stderr,
251  "Points are too concentrated -- please increase DMIN\n");
252  exit(0);
253  }
254 
255  if (data->n_cols % 2 == 0) {
256  cols1 = data->n_cols / 2;
257  cols2 = cols1;
258  }
259  else {
260  cols2 = (int)(data->n_cols / 2);
261  cols1 = cols2 + 1;
262  }
263  if (data->n_rows % 2 == 0) {
264  rows1 = data->n_rows / 2;
265  rows2 = rows1;
266  }
267  else {
268  rows2 = (int)(data->n_rows / 2);
269  rows1 = rows2 + 1;
270  }
271 
272  dx = cols1 * ew_res;
273  dy = rows1 * ns_res;
274 
275  xl = data->x_orig;
276  xm = xl + dx;
277  xr = data->xmax;
278  yl = data->y_orig;
279  ym = yl + dy;
280  yr = data->ymax;
281 
282  if (!(datas = (struct quaddata **)malloc(sizeof(struct quaddata *) * 5))) {
283  return NULL;
284  }
285  datas[NE] = quad_data_new(xm, ym, xr, yr, rows2, cols2, 0, kmax);
286  datas[SW] = quad_data_new(xl, yl, xm, ym, rows1, cols1, 0, kmax);
287  datas[SE] = quad_data_new(xm, yl, xr, ym, rows1, cols2, 0, kmax);
288  datas[NW] = quad_data_new(xl, ym, xm, yr, rows2, cols1, 0, kmax);
289  for (i = 0; i < data->n_points; i++) {
290  switch (quad_compare(data->points + i, data)) {
291  case SW:
292  {
293  quad_add_data(data->points + i, datas[SW], dmin);
294  break;
295  }
296  case SE:
297  {
298  quad_add_data(data->points + i, datas[SE], dmin);
299  break;
300  }
301  case NW:
302  {
303  quad_add_data(data->points + i, datas[NW], dmin);
304  break;
305  }
306  case NE:
307  {
308  quad_add_data(data->points + i, datas[NE], dmin);
309  break;
310  }
311  }
312  }
313  data->points = NULL;
314  return datas;
315 }
316 
317 
318 /*!
319  * Gets such points from *data* that lie within region determined by
320  * *data_inter*. Called by tree function `region_data()`.
321  */
322 int
323 quad_get_points(struct quaddata *data_inter, struct quaddata *data, int MAX)
324 {
325  int i, ind;
326  int n = 0;
327  int l = 0;
328  double xmin, xmax, ymin, ymax;
329  struct triple *point;
330 
331  xmin = data_inter->x_orig;
332  xmax = data_inter->xmax;
333  ymin = data_inter->y_orig;
334  ymax = data_inter->ymax;
335  for (i = 0; i < data->n_points; i++) {
336  point = data->points + i;
337  if (l >= MAX)
338  return MAX + 1;
339  if ((point->x > xmin) && (point->x < xmax)
340  && (point->y > ymin) && (point->y < ymax)) {
341  ind = data_inter->n_points++;
342  data_inter->points[ind].x = point->x;
343  data_inter->points[ind].y = point->y;
344  data_inter->points[ind].z = point->z;
345  data_inter->points[ind].sm = point->sm;
346  l = l + 1;
347 
348  }
349  }
350  n = l;
351  return (n);
352 }
double y_orig
Definition: dataquad.h:51
#define MAX(a, b)
Definition: kdtree.c:26
int quad_intersect(struct quaddata *data_inter, struct quaddata *data)
Definition: dataquad.c:180
double z
Definition: dataquad.h:44
#define NE
Definition: dataquad.h:30
#define NW
Definition: dataquad.h:29
#define NULL
Definition: ccmath.h:32
double x_orig
Definition: dataquad.h:50
#define x
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition: dataquad.c:38
struct triple * points
Definition: dataquad.h:57
void * malloc(YYSIZE_T)
#define SW
Definition: dataquad.h:31
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:62
int quad_compare(struct triple *point, struct quaddata *data)
Definition: dataquad.c:98
double l
Definition: r_raster.c:39
int quad_get_points(struct quaddata *data_inter, struct quaddata *data, int MAX)
Definition: dataquad.c:323
struct quaddata ** quad_divide_data(struct quaddata *data, int kmax, double dmin)
Definition: dataquad.c:237
double x
Definition: dataquad.h:42
double sm
Definition: dataquad.h:45
double ymax
Definition: dataquad.h:53
int quad_division_check(struct quaddata *data, int kmax)
Definition: dataquad.c:219
#define SE
Definition: dataquad.h:32
double xmax
Definition: dataquad.h:52
int n_rows
Definition: dataquad.h:54
int n_cols
Definition: dataquad.h:55
int quad_add_data(struct triple *point, struct quaddata *data, double dmin)
Definition: dataquad.c:143
double y
Definition: dataquad.h:43
double r
Definition: r_raster.c:39
int n_points
Definition: dataquad.h:56