GRASS 8 Programmer's Manual 8.6.0dev(2026)-f6f2c534ea
Loading...
Searching...
No Matches
rect.c
Go to the documentation of this file.
1/****************************************************************************
2 * MODULE: R-Tree library
3 *
4 * AUTHOR(S): Antonin Guttman - original code
5 * Daniel Green (green@superliminal.com) - major clean-up
6 * and implementation of bounding spheres
7 * Markus Metz - file-based and memory-based R*-tree
8 *
9 * PURPOSE: Multidimensional index
10 *
11 * COPYRIGHT: (C) 2010 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *****************************************************************************/
17
18#include <stdio.h>
19#include <stdlib.h>
20#include <assert.h>
21#include "index.h"
22
23#include <float.h>
24#include <math.h>
25#include <grass/gis.h>
26
27#define BIG_NUM (FLT_MAX / 4.0)
28
29#define Undefined(x, t) ((x)->boundary[0] > (x)->boundary[t->ndims_alloc])
30
31/*!
32 \brief Create a new rectangle for a given tree
33
34 This method allocates a new rectangle and initializes
35 the internal boundary coordinates based on the tree dimension.
36
37 Hence a call to RTreeNewBoundary() is not necessary.
38
39 \param t The pointer to a RTree struct
40 \return A new allocated RTree_Rect struct
41 */
43{
44 struct RTree_Rect *r;
45
46 assert(t);
47
48 r = (struct RTree_Rect *)malloc(sizeof(struct RTree_Rect));
49
50 assert(r);
51
52 r->boundary = RTreeAllocBoundary(t);
53 return r;
54}
55
56/*!
57 \brief Delete a rectangle
58
59 This method deletes (free) the allocated memory of a rectangle.
60
61 \param r The pointer to the rectangle to be deleted
62 */
64{
65 assert(r);
67 free(r);
68}
69
70/*!
71 \brief Allocate the boundary array of a rectangle for a given tree
72
73 This method allocated the boundary coordinates array in
74 provided rectangle. It does not release previously allocated memory.
75
76 \param r The pointer to rectangle to initialize the boundary coordinates.
77 This is usually a rectangle that was created on the stack or
78 self allocated.
79 \param t The pointer to a RTree struct
80 */
82{
83 RectReal *boundary = (RectReal *)malloc(t->rectsize);
84
86
87 return boundary;
88}
89
90/*!
91 \brief Delete the boundary of a rectangle
92
93 This method deletes (free) the memory of the boundary of a rectangle
94 and sets the boundary pointer to NULL.
95
96 \param r The pointer to the rectangle to delete the boundary from.
97 */
99{
100 assert(r);
101 if (r->boundary)
102 free(r->boundary);
103 r->boundary = NULL;
104}
105
106/*!
107 \brief Initialize a rectangle to have all 0 coordinates.
108 */
109void RTreeInitRect(struct RTree_Rect *r, struct RTree *t)
110{
111 register int i;
112
113 for (i = 0; i < t->ndims_alloc; i++)
114 r->boundary[i] = r->boundary[i + t->ndims_alloc] = (RectReal)0;
115}
116
117/*!
118 \brief Set one dimensional coordinates of a rectangle for a given tree.
119
120 All coordinates of the rectangle will be initialized to 0 before
121 the x coordinates are set.
122
123 \param r The pointer to the rectangle
124 \param t The pointer to the RTree
125 \param x_min The lower x coordinate
126 \param x_max The higher x coordinate
127 */
128void RTreeSetRect1D(struct RTree_Rect *r, struct RTree *t, double x_min,
129 double x_max)
130{
131 RTreeInitRect(r, t);
132 r->boundary[0] = (RectReal)x_min;
133 r->boundary[t->ndims_alloc] = (RectReal)x_max;
134}
135
136/*!
137 \brief Set two dimensional coordinates of a rectangle for a given tree.
138
139 All coordinates of the rectangle will be initialized to 0 before
140 the x and y coordinates are set.
141
142 \param r The pointer to the rectangle
143 \param t The pointer to the RTree
144 \param x_min The lower x coordinate
145 \param x_max The higher x coordinate
146 \param y_min The lower y coordinate
147 \param y_max The higher y coordinate
148 */
149void RTreeSetRect2D(struct RTree_Rect *r, struct RTree *t, double x_min,
150 double x_max, double y_min, double y_max)
151{
152 RTreeInitRect(r, t);
153 r->boundary[0] = (RectReal)x_min;
154 r->boundary[t->ndims_alloc] = (RectReal)x_max;
155 r->boundary[1] = (RectReal)y_min;
156 r->boundary[1 + t->ndims_alloc] = (RectReal)y_max;
157}
158
159/*!
160 \brief Set three dimensional coordinates of a rectangle for a given tree.
161
162 All coordinates of the rectangle will be initialized to 0 before
163 the x,y and z coordinates are set.
164
165 \param r The pointer to the rectangle
166 \param t The pointer to the RTree
167 \param x_min The lower x coordinate
168 \param x_max The higher x coordinate
169 \param y_min The lower y coordinate
170 \param y_max The higher y coordinate
171 \param z_min The lower z coordinate
172 \param z_max The higher z coordinate
173 */
174void RTreeSetRect3D(struct RTree_Rect *r, struct RTree *t, double x_min,
175 double x_max, double y_min, double y_max, double z_min,
176 double z_max)
177{
178 RTreeInitRect(r, t);
179 r->boundary[0] = (RectReal)x_min;
180 r->boundary[t->ndims_alloc] = (RectReal)x_max;
181 r->boundary[1] = (RectReal)y_min;
182 r->boundary[1 + t->ndims_alloc] = (RectReal)y_max;
183 r->boundary[2] = (RectReal)z_min;
184 r->boundary[2 + t->ndims_alloc] = (RectReal)z_max;
185}
186
187/*!
188 \brief Set 4 dimensional coordinates of a rectangle for a given tree.
189
190 All coordinates of the rectangle will be initialized to 0 before
191 the x,y,z and t coordinates are set.
192
193 \param r The pointer to the rectangle
194 \param t The pointer to the RTree
195 \param x_min The lower x coordinate
196 \param x_max The higher x coordinate
197 \param y_min The lower y coordinate
198 \param y_max The higher y coordinate
199 \param z_min The lower z coordinate
200 \param z_max The higher z coordinate
201 \param t_min The lower t coordinate
202 \param t_max The higher t coordinate
203 */
204void RTreeSetRect4D(struct RTree_Rect *r, struct RTree *t, double x_min,
205 double x_max, double y_min, double y_max, double z_min,
206 double z_max, double t_min, double t_max)
207{
208 assert(t->ndims >= 4);
209
210 RTreeInitRect(r, t);
211 r->boundary[0] = (RectReal)x_min;
212 r->boundary[t->ndims_alloc] = (RectReal)x_max;
213 r->boundary[1] = (RectReal)y_min;
214 r->boundary[1 + t->ndims_alloc] = (RectReal)y_max;
215 r->boundary[2] = (RectReal)z_min;
216 r->boundary[2 + t->ndims_alloc] = (RectReal)z_max;
217 r->boundary[3] = (RectReal)t_min;
218 r->boundary[3 + t->ndims_alloc] = (RectReal)t_max;
219}
220
221/*
222 Return a rect whose first low side is higher than its opposite side -
223 interpreted as an undefined rect.
224 */
225void RTreeNullRect(struct RTree_Rect *r, struct RTree *t)
226{
227 register int i;
228
229 /* assert(r); */
230
231 r->boundary[0] = (RectReal)1;
232 r->boundary[t->nsides_alloc - 1] = (RectReal)-1;
233 for (i = 1; i < t->ndims_alloc; i++)
234 r->boundary[i] = r->boundary[i + t->ndims_alloc] = (RectReal)0;
235
236 return;
237}
238
239#if 0
240
241/*
242 Fills in random coordinates in a rectangle.
243 The low side is guaranteed to be less than the high side.
244 */
245void RTreeRandomRect(struct RTree_Rect *R)
246{
247 register struct RTree_Rect *r = R;
248 register int i;
249 register RectReal width;
250
251 for (i = 0; i < NUMDIMS; i++) {
252 /* width from 1 to 1000 / 4, more small ones
253 */
254 width = drand48() * (1000 / 4) + 1;
255
256 /* sprinkle a given size evenly but so they stay in [0,100]
257 */
258 r->boundary[i] = drand48() * (1000 - width); /* low side */
259 r->boundary[i + NUMDIMS] = r->boundary[i] + width; /* high side */
260 }
261}
262
263
264/*
265 Fill in the boundaries for a random search rectangle.
266 Pass in a pointer to a rect that contains all the data,
267 and a pointer to the rect to be filled in.
268 Generated rect is centered randomly anywhere in the data area,
269 and has size from 0 to the size of the data area in each dimension,
270 i.e. search rect can stick out beyond data area.
271 */
272void RTreeSearchRect(struct RTree_Rect *Search, struct RTree_Rect *Data)
273{
274 register struct RTree_Rect *search = Search, *data = Data;
275 register int i, j;
276 register RectReal size, center;
277
278 assert(search);
279 assert(data);
280
281 for (i = 0; i < NUMDIMS; i++) {
282 j = i + NUMDIMS; /* index for high side boundary */
283 if (data->boundary[i] > -BIG_NUM && data->boundary[j] < BIG_NUM) {
284 size = (drand48() * (data->boundary[j] -
285 data->boundary[i] + 1)) / 2;
286 center = data->boundary[i] + drand48() *
287 (data->boundary[j] - data->boundary[i] + 1);
288 search->boundary[i] = center - size / 2;
289 search->boundary[j] = center + size / 2;
290 }
291 else { /* some open boundary, search entire dimension */
292
293 search->boundary[i] = -BIG_NUM;
294 search->boundary[j] = BIG_NUM;
295 }
296 }
297}
298
299#endif
300
301/*
302 Print out the data for a rectangle.
303 */
304void RTreePrintRect(struct RTree_Rect *R, int depth, struct RTree *t)
305{
306 register struct RTree_Rect *r = R;
307 register int i;
308
309 assert(r);
310
311 RTreeTabIn(depth);
312 fprintf(stdout, "rect:\n");
313 for (i = 0; i < t->ndims_alloc; i++) {
314 RTreeTabIn(depth + 1);
315 fprintf(stdout, "%f\t%f\n", r->boundary[i],
316 r->boundary[i + t->ndims_alloc]);
317 }
318}
319
320/*
321 Calculate the n-dimensional volume of a rectangle
322 */
324{
325 register struct RTree_Rect *r = R;
326 register int i;
327 register RectReal volume = (RectReal)1;
328
329 /* assert(r); */
330
331 if (Undefined(r, t))
332 return (RectReal)0;
333
334 for (i = 0; i < t->ndims; i++)
335 volume *= r->boundary[i + t->ndims_alloc] - r->boundary[i];
336 assert(volume >= 0.0);
337
338 return volume;
339}
340
341/*
342 Define the NUMDIMS-dimensional volume the unit sphere in that dimension into
343 the symbol "UnitSphereVolume"
344 Note that if the gamma function is available in the math library and if the
345 compiler supports static initialization using functions, this is
346 easily computed for any dimension. If not, the value can be precomputed and
347 taken from a table. The following code can do it either way.
348 */
349
350#ifdef gamma
351
352/* computes the volume of an N-dimensional sphere. */
353/* derived from formula in "Regular Polytopes" by H.S.M Coxeter */
354static double sphere_volume(double dimension)
355{
356 double log_gamma, log_volume;
357
358 log_gamma = gamma(dimension / 2.0 + 1);
359 log_volume = dimension / 2.0 * log(M_PI) - log_gamma;
360 return exp(log_volume);
361}
362
363static const double UnitSphereVolume = sphere_volume(20);
364
365#else
366
367/* Precomputed volumes of the unit spheres for the first few dimensions */
368const double UnitSphereVolumes[] = {
369 0.000000, /* dimension 0 */
370 2.000000, /* dimension 1 */
371 3.141593, /* dimension 2 */
372 4.188790, /* dimension 3 */
373 4.934802, /* dimension 4 */
374 5.263789, /* dimension 5 */
375 5.167713, /* dimension 6 */
376 4.724766, /* dimension 7 */
377 4.058712, /* dimension 8 */
378 3.298509, /* dimension 9 */
379 2.550164, /* dimension 10 */
380 1.884104, /* dimension 11 */
381 1.335263, /* dimension 12 */
382 0.910629, /* dimension 13 */
383 0.599265, /* dimension 14 */
384 0.381443, /* dimension 15 */
385 0.235331, /* dimension 16 */
386 0.140981, /* dimension 17 */
387 0.082146, /* dimension 18 */
388 0.046622, /* dimension 19 */
389 0.025807, /* dimension 20 */
390};
391
392#if NUMDIMS > 20
393#error "not enough precomputed sphere volumes"
394#endif
395#define UnitSphereVolume UnitSphereVolumes[NUMDIMS]
396
397#endif
398
399/*
400 Calculate the n-dimensional volume of the bounding sphere of a rectangle
401 */
402
403#if 0
404/*
405 * A fast approximation to the volume of the bounding sphere for the
406 * given Rect. By Paul B.
407 */
409{
410 register struct RTree_Rect *r = R;
411 register int i;
412 RectReal maxsize = (RectReal) 0, c_size;
413
414 /* assert(r); */
415
416 if (Undefined(r, t))
417 return (RectReal) 0;
418
419 for (i = 0; i < t->ndims; i++) {
420 c_size = r->boundary[i + NUMDIMS] - r->boundary[i];
421 if (c_size > maxsize)
422 maxsize = c_size;
423 }
424 return (RectReal) (pow(maxsize / 2, NUMDIMS) *
425 UnitSphereVolumes[t->ndims]);
426}
427#endif
428
429/*
430 * The exact volume of the bounding sphere for the given Rect.
431 */
433{
434 int i;
435 double sum_of_squares = 0, extent;
436
437 /* assert(r); */
438
439 if (Undefined(r, t))
440 return (RectReal)0;
441
442 for (i = 0; i < t->ndims; i++) {
443 extent = (r->boundary[i + t->ndims_alloc] - r->boundary[i]);
444
445 /* extent should be half extent : /4 */
446 sum_of_squares += extent * extent / 4.;
447 }
448
449 return (RectReal)(pow(sqrt(sum_of_squares), t->ndims) *
450 UnitSphereVolumes[t->ndims]);
451}
452
453/*
454 Calculate the n-dimensional surface area of a rectangle
455 */
457{
458 int i, j;
459 RectReal face_area, sum = (RectReal)0;
460
461 /*assert(r); */
462
463 if (Undefined(r, t))
464 return (RectReal)0;
465
466 for (i = 0; i < t->ndims; i++) {
467 face_area = (RectReal)1;
468
469 for (j = 0; j < t->ndims; j++)
470 /* exclude i extent from product in this dimension */
471 if (i != j) {
472 face_area *= (r->boundary[j + t->ndims_alloc] - r->boundary[j]);
473 }
474 sum += face_area;
475 }
476 return 2 * sum;
477}
478
479/*
480 Calculate the n-dimensional margin of a rectangle
481 the margin is the sum of the lengths of the edges
482 */
484{
485 int i;
486 RectReal margin = 0.0;
487
488 /* assert(r); */
489
490 for (i = 0; i < t->ndims; i++) {
491 margin += r->boundary[i + t->ndims_alloc] - r->boundary[i];
492 }
493
494 return margin;
495}
496
497/*
498 Combine two rectangles, make one that includes both.
499 */
501 struct RTree_Rect *r3, struct RTree *t)
502{
503 int i, j;
504
505 /* assert(r1 && r2 && r3); */
506
507 if (Undefined(r1, t)) {
508 for (i = 0; i < t->nsides_alloc; i++)
509 r3->boundary[i] = r2->boundary[i];
510
511 return;
512 }
513
514 if (Undefined(r2, t)) {
515 for (i = 0; i < t->nsides_alloc; i++)
516 r3->boundary[i] = r1->boundary[i];
517
518 return;
519 }
520
521 for (i = 0; i < t->ndims; i++) {
522 r3->boundary[i] = MIN(r1->boundary[i], r2->boundary[i]);
523 j = i + t->ndims_alloc;
524 r3->boundary[j] = MAX(r1->boundary[j], r2->boundary[j]);
525 }
526 for (i = t->ndims; i < t->ndims_alloc; i++) {
527 r3->boundary[i] = 0;
528 j = i + t->ndims_alloc;
529 r3->boundary[j] = 0;
530 }
531}
532
533/*
534 Expand first rectangle to cover second rectangle.
535 */
537 struct RTree *t)
538{
539 int i, j, ret = 0;
540
541 /* assert(r1 && r2); */
542
543 if (Undefined(r2, t))
544 return ret;
545
546 for (i = 0; i < t->ndims; i++) {
547 if (r1->boundary[i] > r2->boundary[i]) {
548 r1->boundary[i] = r2->boundary[i];
549 ret = 1;
550 }
551 j = i + t->ndims_alloc;
552 if (r1->boundary[j] < r2->boundary[j]) {
553 r1->boundary[j] = r2->boundary[j];
554 ret = 1;
555 }
556 }
557
558 for (i = t->ndims; i < t->ndims_alloc; i++) {
559 r1->boundary[i] = 0;
560 j = i + t->ndims_alloc;
561 r1->boundary[j] = 0;
562 }
563
564 return ret;
565}
566
567/*
568 Decide whether two rectangles are identical.
569 */
571 struct RTree *t)
572{
573 register int i, j;
574
575 /* assert(r && s); */
576
577 for (i = 0; i < t->ndims; i++) {
578 j = i + t->ndims_alloc; /* index for high sides */
579 if (r->boundary[i] != s->boundary[i] ||
580 r->boundary[j] != s->boundary[j]) {
581 return 0;
582 }
583 }
584 return 1;
585}
586
587/*
588 Decide whether two rectangles overlap or touch.
589 */
590int RTreeOverlap(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
591{
592 register int i, j;
593
594 /* assert(r && s); */
595
596 for (i = 0; i < t->ndims; i++) {
597 j = i + t->ndims_alloc; /* index for high sides */
598 if (r->boundary[i] > s->boundary[j] ||
599 s->boundary[i] > r->boundary[j]) {
600 return FALSE;
601 }
602 }
603 return TRUE;
604}
605
606/*
607 Decide whether rectangle s is contained in rectangle r.
608 */
609int RTreeContained(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
610{
611 register int i, j;
612
613 /* assert(r && s); */
614
615 /* undefined rect is contained in any other */
616 if (Undefined(r, t))
617 return TRUE;
618
619 /* no rect (except an undefined one) is contained in an undef rect */
620 if (Undefined(s, t))
621 return FALSE;
622
623 for (i = 0; i < t->ndims; i++) {
624 j = i + t->ndims_alloc; /* index for high sides */
625 if (s->boundary[i] < r->boundary[i] || s->boundary[j] > r->boundary[j])
626 return FALSE;
627 }
628 return TRUE;
629}
630
631/*
632 Decide whether rectangle s fully contains rectangle r.
633 */
634int RTreeContains(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
635{
636 register int i, j;
637
638 /* assert(r && s); */
639
640 /* undefined rect is contained in any other */
641 if (Undefined(r, t))
642 return TRUE;
643
644 /* no rect (except an undefined one) is contained in an undef rect */
645 if (Undefined(s, t))
646 return FALSE;
647
648 for (i = 0; i < t->ndims; i++) {
649 j = i + t->ndims_alloc; /* index for high sides */
650 if (s->boundary[i] > r->boundary[i] || s->boundary[j] < r->boundary[j])
651 return FALSE;
652 }
653 return TRUE;
654}
#define NULL
Definition ccmath.h:32
double sphere_volume(double dimension)
Definition gammavol.c:27
#define MIN(a, b)
Definition gis.h:153
#define TRUE
Definition gis.h:78
#define FALSE
Definition gis.h:82
#define M_PI
Definition gis.h:157
#define MAX(a, b)
Definition gis.h:148
void RTreeTabIn(int)
Definition node.c:602
#define assert(condition)
Definition lz4.c:291
double t
Definition r_raster.c:39
double r
Definition r_raster.c:39
struct RTree_Rect * RTreeAllocRect(struct RTree *t)
Create a new rectangle for a given tree.
Definition rect.c:42
void RTreeFreeRect(struct RTree_Rect *r)
Delete a rectangle.
Definition rect.c:63
void RTreeNullRect(struct RTree_Rect *r, struct RTree *t)
Definition rect.c:225
const double UnitSphereVolumes[]
Definition rect.c:368
#define Undefined(x, t)
Definition rect.c:29
void RTreeCombineRect(struct RTree_Rect *r1, struct RTree_Rect *r2, struct RTree_Rect *r3, struct RTree *t)
Definition rect.c:500
#define BIG_NUM
Definition rect.c:27
RectReal RTreeRectVolume(struct RTree_Rect *R, struct RTree *t)
Definition rect.c:323
void RTreeFreeBoundary(struct RTree_Rect *r)
Delete the boundary of a rectangle.
Definition rect.c:98
void RTreeSetRect1D(struct RTree_Rect *r, struct RTree *t, double x_min, double x_max)
Set one dimensional coordinates of a rectangle for a given tree.
Definition rect.c:128
void RTreeSetRect2D(struct RTree_Rect *r, struct RTree *t, double x_min, double x_max, double y_min, double y_max)
Set two dimensional coordinates of a rectangle for a given tree.
Definition rect.c:149
void RTreeSetRect4D(struct RTree_Rect *r, struct RTree *t, double x_min, double x_max, double y_min, double y_max, double z_min, double z_max, double t_min, double t_max)
Set 4 dimensional coordinates of a rectangle for a given tree.
Definition rect.c:204
void RTreeInitRect(struct RTree_Rect *r, struct RTree *t)
Initialize a rectangle to have all 0 coordinates.
Definition rect.c:109
RectReal * RTreeAllocBoundary(struct RTree *t)
Allocate the boundary array of a rectangle for a given tree.
Definition rect.c:81
int RTreeContains(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
Definition rect.c:634
RectReal RTreeRectSurfaceArea(struct RTree_Rect *r, struct RTree *t)
Definition rect.c:456
int RTreeExpandRect(struct RTree_Rect *r1, struct RTree_Rect *r2, struct RTree *t)
Definition rect.c:536
RectReal RTreeRectMargin(struct RTree_Rect *r, struct RTree *t)
Definition rect.c:483
int RTreeOverlap(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
Definition rect.c:590
int RTreeContained(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
Definition rect.c:609
int RTreeCompareRect(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
Definition rect.c:570
#define UnitSphereVolume
Definition rect.c:395
void RTreePrintRect(struct RTree_Rect *R, int depth, struct RTree *t)
Definition rect.c:304
RectReal RTreeRectSphericalVolume(struct RTree_Rect *r, struct RTree *t)
Definition rect.c:432
void RTreeSetRect3D(struct RTree_Rect *r, struct RTree *t, double x_min, double x_max, double y_min, double y_max, double z_min, double z_max)
Set three dimensional coordinates of a rectangle for a given tree.
Definition rect.c:174
double RectReal
Definition rtree.h:26
void * malloc(unsigned)
void free(void *)
RectReal * boundary
Definition rtree.h:55
Definition rtree.h:123