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