GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
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 */
38struct 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 */
62struct 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 free(data);
84 return NULL;
85 }
86 for (i = 0; i <= kmax; i++) {
87 data->points[i].x = 0.;
88 data->points[i].y = 0.;
89 data->points[i].z = 0.;
90 data->points[i].sm = 0.;
91 }
92
93 return data;
94}
95
96
97/*!
98 * Return the quadrant the point should be inserted in
99 */
100int quad_compare(struct triple *point, struct quaddata *data)
101{
102 int cond1, cond2, cond3, cond4, rows, cols;
103 double ew_res, ns_res;
104
105 if (data == NULL)
106 return -1;
107
108 ew_res = (data->xmax - data->x_orig) / data->n_cols;
109 ns_res = (data->ymax - data->y_orig) / data->n_rows;
110
111 if (data->n_rows % 2 == 0) {
112 rows = data->n_rows / 2;
113 }
114 else {
115 rows = (int)(data->n_rows / 2) + 1;
116 }
117
118 if (data->n_cols % 2 == 0) {
119 cols = data->n_cols / 2;
120 }
121 else {
122 cols = (int)(data->n_cols / 2) + 1;
123 }
124 cond1 = (point->x >= data->x_orig);
125 cond2 = (point->x >= data->x_orig + ew_res * cols);
126 cond3 = (point->y >= data->y_orig);
127 cond4 = (point->y >= data->y_orig + ns_res * rows);
128 if (cond1 && cond3) {
129 if (cond2 && cond4)
130 return NE;
131 if (cond2)
132 return SE;
133 if (cond4)
134 return NW;
135 return SW;
136 }
137 else
138 return 0;
139}
140
141
142/*!
143 * Add point to a given *data*.
144 */
145int quad_add_data(struct triple *point, struct quaddata *data, double dmin)
146{
147
148 int cond = 1;
149 if (data == NULL) {
150 fprintf(stderr, "add_data: data is NULL \n");
151 return -5;
152 }
153 for (int i = 0; i < data->n_points; i++) {
154 double xx = data->points[i].x - point->x;
155 double yy = data->points[i].y - point->y;
156 double r = xx * xx + yy * yy;
157 if (r <= dmin) {
158 cond = 0;
159 break;
160 }
161 }
162
163 if (cond) {
164 int 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 */
180int 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 */
219int 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 */
237struct 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 */
322int
323quad_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}
#define NULL
Definition: ccmath.h:32
int quad_add_data(struct triple *point, struct quaddata *data, double dmin)
Definition: dataquad.c:145
int quad_division_check(struct quaddata *data, int kmax)
Definition: dataquad.c:219
int quad_get_points(struct quaddata *data_inter, struct quaddata *data, int MAX)
Definition: dataquad.c:323
int quad_intersect(struct quaddata *data_inter, struct quaddata *data)
Definition: dataquad.c:180
int quad_compare(struct triple *point, struct quaddata *data)
Definition: dataquad.c:100
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition: dataquad.c:38
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
struct quaddata ** quad_divide_data(struct quaddata *data, int kmax, double dmin)
Definition: dataquad.c:237
#define SE
Definition: dataquad.h:32
#define SW
Definition: dataquad.h:31
#define NE
Definition: dataquad.h:30
#define NW
Definition: dataquad.h:29
double l
double r
#define MAX(a, b)
Definition: shpopen.c:299
double ymax
Definition: dataquad.h:53
double y_orig
Definition: dataquad.h:51
double x_orig
Definition: dataquad.h:50
struct triple * points
Definition: dataquad.h:57
int n_points
Definition: dataquad.h:56
double xmax
Definition: dataquad.h:52
int n_cols
Definition: dataquad.h:55
int n_rows
Definition: dataquad.h:54
double z
Definition: dataquad.h:44
double sm
Definition: dataquad.h:45
double x
Definition: dataquad.h:42
double y
Definition: dataquad.h:43
#define x