GRASS Programmer's Manual  6.5.svn(2014)-r66266
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 *
9 * PURPOSE: Multidimensional index
10 *
11 * COPYRIGHT: (C) 2001 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
30 #define Undefined(x) ((x)->boundary[0] > (x)->boundary[NUMDIMS])
31 #define MIN(a, b) ((a) < (b) ? (a) : (b))
32 #define MAX(a, b) ((a) > (b) ? (a) : (b))
33
34
35 /*-----------------------------------------------------------------------------
36 | Initialize a rectangle to have all 0 coordinates.
37 -----------------------------------------------------------------------------*/
38 void RTreeInitRect(struct Rect *R)
39 {
40  register struct Rect *r = R;
41  register int i;
42
43  for (i = 0; i < NUMSIDES; i++)
44  r->boundary[i] = (RectReal) 0;
45 }
46
47
48 /*-----------------------------------------------------------------------------
49 | Return a rect whose first low side is higher than its opposite side -
50 | interpreted as an undefined rect.
51 -----------------------------------------------------------------------------*/
52 struct Rect RTreeNullRect(void)
53 {
54  struct Rect r;
55  register int i;
56
57  r.boundary[0] = (RectReal) 1;
58  r.boundary[NUMDIMS] = (RectReal) - 1;
59  for (i = 1; i < NUMDIMS; i++)
60  r.boundary[i] = r.boundary[i + NUMDIMS] = (RectReal) 0;
61  return r;
62 }
63
64
65 #if 0
66
67 /*-----------------------------------------------------------------------------
68 | Fills in random coordinates in a rectangle.
69 | The low side is guaranteed to be less than the high side.
70 -----------------------------------------------------------------------------*/
71 void RTreeRandomRect(struct Rect *R)
72 {
73  register struct Rect *r = R;
74  register int i;
75  register RectReal width;
76
77  for (i = 0; i < NUMDIMS; i++) {
78  /* width from 1 to 1000 / 4, more small ones
79  */
80  width = drand48() * (1000 / 4) + 1;
81
82  /* sprinkle a given size evenly but so they stay in [0,100]
83  */
84  r->boundary[i] = drand48() * (1000 - width); /* low side */
85  r->boundary[i + NUMDIMS] = r->boundary[i] + width; /* high side */
86  }
87 }
88
89
90 /*-----------------------------------------------------------------------------
91 | Fill in the boundaries for a random search rectangle.
92 | Pass in a pointer to a rect that contains all the data,
93 | and a pointer to the rect to be filled in.
94 | Generated rect is centered randomly anywhere in the data area,
95 | and has size from 0 to the size of the data area in each dimension,
96 | i.e. search rect can stick out beyond data area.
97 -----------------------------------------------------------------------------*/
98 void RTreeSearchRect(struct Rect *Search, struct Rect *Data)
99 {
100  register struct Rect *search = Search, *data = Data;
101  register int i, j;
102  register RectReal size, center;
103
104  assert(search);
105  assert(data);
106
107  for (i = 0; i < NUMDIMS; i++) {
108  j = i + NUMDIMS; /* index for high side boundary */
109  if (data->boundary[i] > -BIG_NUM && data->boundary[j] < BIG_NUM) {
110  size = (drand48() * (data->boundary[j] -
111  data->boundary[i] + 1)) / 2;
112  center = data->boundary[i] + drand48() *
113  (data->boundary[j] - data->boundary[i] + 1);
114  search->boundary[i] = center - size / 2;
115  search->boundary[j] = center + size / 2;
116  }
117  else { /* some open boundary, search entire dimension */
118
119  search->boundary[i] = -BIG_NUM;
120  search->boundary[j] = BIG_NUM;
121  }
122  }
123 }
124
125 #endif
126
127 /*-----------------------------------------------------------------------------
128 | Print out the data for a rectangle.
129 -----------------------------------------------------------------------------*/
130 void RTreePrintRect(struct Rect *R, int depth)
131 {
132  register struct Rect *r = R;
133  register int i;
134
135  assert(r);
136
137  RTreeTabIn(depth);
138  fprintf(stdout, "rect:\n");
139  for (i = 0; i < NUMDIMS; i++) {
140  RTreeTabIn(depth + 1);
141  fprintf(stdout, "%f\t%f\n", r->boundary[i], r->boundary[i + NUMDIMS]);
142  }
143 }
144
145 /*-----------------------------------------------------------------------------
146 | Calculate the n-dimensional volume of a rectangle
147 -----------------------------------------------------------------------------*/
149 {
150  register struct Rect *r = R;
151  register int i;
152  register RectReal volume = (RectReal) 1;
153
154  assert(r);
155  if (Undefined(r))
156  return (RectReal) 0;
157
158  for (i = 0; i < NUMDIMS; i++)
159  volume *= r->boundary[i + NUMDIMS] - r->boundary[i];
160  assert(volume >= 0.0);
161  return volume;
162 }
163
164
165 /*-----------------------------------------------------------------------------
166 | Define the NUMDIMS-dimensional volume the unit sphere in that dimension into
167 | the symbol "UnitSphereVolume"
168 | Note that if the gamma function is available in the math library and if the
169 | compiler supports static initialization using functions, this is
170 | easily computed for any dimension. If not, the value can be precomputed and
171 | taken from a table. The following code can do it either way.
172 -----------------------------------------------------------------------------*/
173
174 #ifdef gamma
175
176 /* computes the volume of an N-dimensional sphere. */
177 /* derived from formule in "Regular Polytopes" by H.S.M Coxeter */
178 static double sphere_volume(double dimension)
179 {
180  double log_gamma, log_volume;
181
182  log_gamma = gamma(dimension / 2.0 + 1);
183  log_volume = dimension / 2.0 * log(M_PI) - log_gamma;
184  return exp(log_volume);
185 }
186 static const double UnitSphereVolume = sphere_volume(NUMDIMS);
187
188 #else
189
190 /* Precomputed volumes of the unit spheres for the first few dimensions */
191 const double UnitSphereVolumes[] = {
192  0.000000, /* dimension 0 */
193  2.000000, /* dimension 1 */
194  3.141593, /* dimension 2 */
195  4.188790, /* dimension 3 */
196  4.934802, /* dimension 4 */
197  5.263789, /* dimension 5 */
198  5.167713, /* dimension 6 */
199  4.724766, /* dimension 7 */
200  4.058712, /* dimension 8 */
201  3.298509, /* dimension 9 */
202  2.550164, /* dimension 10 */
203  1.884104, /* dimension 11 */
204  1.335263, /* dimension 12 */
205  0.910629, /* dimension 13 */
206  0.599265, /* dimension 14 */
207  0.381443, /* dimension 15 */
208  0.235331, /* dimension 16 */
209  0.140981, /* dimension 17 */
210  0.082146, /* dimension 18 */
211  0.046622, /* dimension 19 */
212  0.025807, /* dimension 20 */
213 };
214
215 #if NUMDIMS > 20
216 # error "not enough precomputed sphere volumes"
217 #endif
218 #define UnitSphereVolume UnitSphereVolumes[NUMDIMS]
219
220 #endif
221
222
223 /*-----------------------------------------------------------------------------
224 | Calculate the n-dimensional volume of the bounding sphere of a rectangle
225 -----------------------------------------------------------------------------*/
226
227 #if 0
228 /*
229  * A fast approximation to the volume of the bounding sphere for the
230  * given Rect. By Paul B.
231  */
233 {
234  register struct Rect *r = R;
235  register int i;
236  RectReal maxsize = (RectReal) 0, c_size;
237
238  assert(r);
239  if (Undefined(r))
240  return (RectReal) 0;
241  for (i = 0; i < NUMDIMS; i++) {
242  c_size = r->boundary[i + NUMDIMS] - r->boundary[i];
243  if (c_size > maxsize)
244  maxsize = c_size;
245  }
246  return (RectReal) (pow(maxsize / 2, NUMDIMS) * UnitSphereVolume);
247 }
248 #endif
249
250 /*
251  * The exact volume of the bounding sphere for the given Rect.
252  */
254 {
255  register struct Rect *r = R;
256  register int i;
257  register double sum_of_squares = 0, radius;
258
259  assert(r);
260  if (Undefined(r))
261  return (RectReal) 0;
262  for (i = 0; i < NUMDIMS; i++) {
263  double half_extent = (r->boundary[i + NUMDIMS] - r->boundary[i]) / 2;
264
265  sum_of_squares += half_extent * half_extent;
266  }
268  return (RectReal) (pow(radius, NUMDIMS) * UnitSphereVolume);
269 }
270
271
272 /*-----------------------------------------------------------------------------
273 | Calculate the n-dimensional surface area of a rectangle
274 -----------------------------------------------------------------------------*/
276 {
277  register struct Rect *r = R;
278  register int i, j;
279  register RectReal sum = (RectReal) 0;
280
281  assert(r);
282  if (Undefined(r))
283  return (RectReal) 0;
284
285  for (i = 0; i < NUMDIMS; i++) {
286  RectReal face_area = (RectReal) 1;
287
288  for (j = 0; j < NUMDIMS; j++)
289  /* exclude i extent from product in this dimension */
290  if (i != j) {
291  RectReal j_extent = r->boundary[j + NUMDIMS] - r->boundary[j];
292
293  face_area *= j_extent;
294  }
295  sum += face_area;
296  }
297  return 2 * sum;
298 }
299
300
301
302 /*-----------------------------------------------------------------------------
303 | Combine two rectangles, make one that includes both.
304 -----------------------------------------------------------------------------*/
305 struct Rect RTreeCombineRect(struct Rect *R, struct Rect *Rr)
306 {
307  register struct Rect *r = R, *rr = Rr;
308  register int i, j;
309  struct Rect new_rect;
310
311  assert(r && rr);
312
313  if (Undefined(r))
314  return *rr;
315
316  if (Undefined(rr))
317  return *r;
318
319  for (i = 0; i < NUMDIMS; i++) {
320  new_rect.boundary[i] = MIN(r->boundary[i], rr->boundary[i]);
321  j = i + NUMDIMS;
322  new_rect.boundary[j] = MAX(r->boundary[j], rr->boundary[j]);
323  }
324  return new_rect;
325 }
326
327
328 /*-----------------------------------------------------------------------------
329 | Decide whether two rectangles overlap.
330 -----------------------------------------------------------------------------*/
331 int RTreeOverlap(struct Rect *R, struct Rect *S)
332 {
333  register struct Rect *r = R, *s = S;
334  register int i, j;
335
336  assert(r && s);
337
338  for (i = 0; i < NUMDIMS; i++) {
339  j = i + NUMDIMS; /* index for high sides */
340  if (r->boundary[i] > s->boundary[j] ||
341  s->boundary[i] > r->boundary[j]) {
342  return FALSE;
343  }
344  }
345  return TRUE;
346 }
347
348
349 /*-----------------------------------------------------------------------------
350 | Decide whether rectangle r is contained in rectangle s.
351 -----------------------------------------------------------------------------*/
352 int RTreeContained(struct Rect *R, struct Rect *S)
353 {
354  register struct Rect *r = R, *s = S;
355  register int i, j, result;
356
357  assert(r && s); /* same as in RTreeOverlap() */
358
359  /* undefined rect is contained in any other */
360  if (Undefined(r))
361  return TRUE;
362
363  /* no rect (except an undefined one) is contained in an undef rect */
364  if (Undefined(s))
365  return FALSE;
366
367  result = TRUE;
368  for (i = 0; i < NUMDIMS; i++) {
369  j = i + NUMDIMS; /* index for high sides */
370  result = result && r->boundary[i] >= s->boundary[i]
371  && r->boundary[j] <= s->boundary[j];
372  }
373  return result;
374 }
RectReal boundary[NUMSIDES]
Definition: index.h:42
double RectReal
Definition: index.h:25
#define FALSE
Definition: dbfopen.c:117
#define UnitSphereVolume
Definition: rect.c:218
#define NUMSIDES
Definition: index.h:38
#define MAX(a, b)
Definition: rect.c:32
float r
Definition: named_colr.c:8
tuple width
RectReal RTreeRectSphericalVolume(struct Rect *R)
Definition: rect.c:253
int RTreeContained(struct Rect *R, struct Rect *S)
Definition: rect.c:352
#define MIN(a, b)
Definition: rect.c:31
#define BIG_NUM
Definition: rect.c:27
#define NUMDIMS
Definition: index.h:22
RectReal RTreeRectVolume(struct Rect *R)
Definition: rect.c:148
tuple size
value.Bind(wx.EVT_TEXT, self.OnVolumeIsosurfMap)
Definition: tools.py:2334
double sphere_volume(double dimension)
Definition: gammavol.c:27
log
Definition: wxnviz.py:54
tuple data
void RTreeInitRect(struct Rect *)
Definition: rect.c:38
int RTreeOverlap(struct Rect *, struct Rect *)
Definition: rect.c:331
#define TRUE
Definition: dbfopen.c:118
void RTreePrintRect(struct Rect *, int)
Definition: rect.c:130
const double UnitSphereVolumes[]
Definition: rect.c:191
struct Rect RTreeNullRect(void)
Definition: rect.c:52
#define Undefined(x)
Definition: rect.c:30
struct Rect RTreeCombineRect(struct Rect *, struct Rect *)
Definition: rect.c:305
RectReal RTreeRectSurfaceArea(struct Rect *R)
Definition: rect.c:275
Definition: index.h:40
void RTreeTabIn(int)
Definition: node.c:68