GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-c1f55e303d
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  */
42 struct RTree_Rect *RTreeAllocRect(struct RTree *t)
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  */
63 void RTreeFreeRect(struct RTree_Rect *r)
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  */
109 void 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  */
128 void 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  */
149 void 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  */
174 void 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  */
204 void 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  */
225 void 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  */
245 void 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  */
272 void 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  */
304 void 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 formule in "Regular Polytopes" by H.S.M Coxeter */
354 static 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
363 static const double UnitSphereVolume = sphere_volume(20);
364
365 #else
366
367 /* Precomputed volumes of the unit spheres for the first few dimensions */
368 const 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  */
500 void RTreeCombineRect(struct RTree_Rect *r1, struct RTree_Rect *r2,
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  */
536 int RTreeExpandRect(struct RTree_Rect *r1, struct RTree_Rect *r2,
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  */
570 int RTreeCompareRect(struct RTree_Rect *r, struct RTree_Rect *s,
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  */
590 int 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  */
609 int 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  */
634 int 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:154
#define TRUE
Definition: gis.h:79
#define FALSE
Definition: gis.h:83
#define M_PI
Definition: gis.h:158
#define MAX(a, b)
Definition: gis.h:149
void RTreeTabIn(int)
Definition: node.c:604
#define assert(condition)
Definition: lz4.c:393
double t
Definition: r_raster.c:39
double r
Definition: r_raster.c:39
RectReal * RTreeAllocBoundary(struct RTree *t)
Allocate the boundary array of a rectangle for a given tree.
Definition: rect.c:81
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
struct RTree_Rect * RTreeAllocRect(struct RTree *t)
Create a new rectangle for a given tree.
Definition: rect.c:42
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
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(YYSIZE_T)
void free(void *)
RectReal * boundary
Definition: rtree.h:55
Definition: rtree.h:123