GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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 */
36struct 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 */
59struct 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 */
95int 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 */
139int 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 */
176int 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 */
209int 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 */
226struct 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 */
305int 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
int quad_compare(struct triple *point, struct quaddata *data)
Definition dataquad.c:95
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition dataquad.c:36
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
struct quaddata ** quad_divide_data(struct quaddata *data, int kmax, double dmin)
Definition dataquad.c:226
#define SE
Definition dataquad.h:31
#define SW
Definition dataquad.h:30
#define NE
Definition dataquad.h:29
#define NW
Definition dataquad.h:28
double l
double r
#define MAX(a, b)
Definition shpopen.c:65
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