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