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