/**************************************************************************** * MODULE: R-Tree library * * AUTHOR(S): Antonin Guttman - original code * Daniel Green (green@superliminal.com) - major clean-up * and implementation of bounding spheres * Markus Metz - R*-tree * * PURPOSE: Multidimensional index * * COPYRIGHT: (C) 2009 by the GRASS Development Team * * This program is free software under the GNU General Public * License (>=v2). Read the file COPYING that comes with GRASS * for details. *****************************************************************************/ #include #include #include #include "index.h" #include #include #define BIG_NUM (FLT_MAX/4.0) #define Undefined(x) ((x)->boundary[0] > (x)->boundary[NUMDIMS]) #define MIN(a, b) ((a) < (b) ? (a) : (b)) #define MAX(a, b) ((a) > (b) ? (a) : (b)) /*----------------------------------------------------------------------------- | Initialize a rectangle to have all 0 coordinates. -----------------------------------------------------------------------------*/ void RTreeInitRect(struct Rect *R) { register struct Rect *r = R; register int i; for (i = 0; i < NUMSIDES; i++) r->boundary[i] = (RectReal) 0; } /*----------------------------------------------------------------------------- | Return a rect whose first low side is higher than its opposite side - | interpreted as an undefined rect. -----------------------------------------------------------------------------*/ struct Rect RTreeNullRect(void) { struct Rect r; register int i; r.boundary[0] = (RectReal) 1; r.boundary[NUMDIMS] = (RectReal) - 1; for (i = 1; i < NUMDIMS; i++) r.boundary[i] = r.boundary[i + NUMDIMS] = (RectReal) 0; return r; } #if 0 /*----------------------------------------------------------------------------- | Fills in random coordinates in a rectangle. | The low side is guaranteed to be less than the high side. -----------------------------------------------------------------------------*/ void RTreeRandomRect(struct Rect *R) { register struct Rect *r = R; register int i; register RectReal width; for (i = 0; i < NUMDIMS; i++) { /* width from 1 to 1000 / 4, more small ones */ width = drand48() * (1000 / 4) + 1; /* sprinkle a given size evenly but so they stay in [0,100] */ r->boundary[i] = drand48() * (1000 - width); /* low side */ r->boundary[i + NUMDIMS] = r->boundary[i] + width; /* high side */ } } /*----------------------------------------------------------------------------- | Fill in the boundaries for a random search rectangle. | Pass in a pointer to a rect that contains all the data, | and a pointer to the rect to be filled in. | Generated rect is centered randomly anywhere in the data area, | and has size from 0 to the size of the data area in each dimension, | i.e. search rect can stick out beyond data area. -----------------------------------------------------------------------------*/ void RTreeSearchRect(struct Rect *Search, struct Rect *Data) { register struct Rect *search = Search, *data = Data; register int i, j; register RectReal size, center; assert(search); assert(data); for (i = 0; i < NUMDIMS; i++) { j = i + NUMDIMS; /* index for high side boundary */ if (data->boundary[i] > -BIG_NUM && data->boundary[j] < BIG_NUM) { size = (drand48() * (data->boundary[j] - data->boundary[i] + 1)) / 2; center = data->boundary[i] + drand48() * (data->boundary[j] - data->boundary[i] + 1); search->boundary[i] = center - size / 2; search->boundary[j] = center + size / 2; } else { /* some open boundary, search entire dimension */ search->boundary[i] = -BIG_NUM; search->boundary[j] = BIG_NUM; } } } #endif /*----------------------------------------------------------------------------- | Print out the data for a rectangle. -----------------------------------------------------------------------------*/ void RTreePrintRect(struct Rect *R, int depth) { register struct Rect *r = R; register int i; assert(r); RTreeTabIn(depth); fprintf(stdout, "rect:\n"); for (i = 0; i < NUMDIMS; i++) { RTreeTabIn(depth + 1); fprintf(stdout, "%f\t%f\n", r->boundary[i], r->boundary[i + NUMDIMS]); } } /*----------------------------------------------------------------------------- | Calculate the n-dimensional volume of a rectangle -----------------------------------------------------------------------------*/ RectReal RTreeRectVolume(struct Rect *R, struct RTree *t) { register struct Rect *r = R; register int i; register RectReal volume = (RectReal) 1; assert(r); if (Undefined(r)) return (RectReal) 0; for (i = 0; i < t->ndims; i++) volume *= r->boundary[i + NUMDIMS] - r->boundary[i]; assert(volume >= 0.0); return volume; } /*----------------------------------------------------------------------------- | Define the NUMDIMS-dimensional volume the unit sphere in that dimension into | the symbol "UnitSphereVolume" | Note that if the gamma function is available in the math library and if the | compiler supports static initialization using functions, this is | easily computed for any dimension. If not, the value can be precomputed and | taken from a table. The following code can do it either way. -----------------------------------------------------------------------------*/ #ifdef gamma /* computes the volume of an N-dimensional sphere. */ /* derived from formule in "Regular Polytopes" by H.S.M Coxeter */ static double sphere_volume(double dimension) { double log_gamma, log_volume; log_gamma = gamma(dimension / 2.0 + 1); log_volume = dimension / 2.0 * log(M_PI) - log_gamma; return exp(log_volume); } static const double UnitSphereVolume = sphere_volume(NUMDIMS); #else /* Precomputed volumes of the unit spheres for the first few dimensions */ const double UnitSphereVolumes[] = { 0.000000, /* dimension 0 */ 2.000000, /* dimension 1 */ 3.141593, /* dimension 2 */ 4.188790, /* dimension 3 */ 4.934802, /* dimension 4 */ 5.263789, /* dimension 5 */ 5.167713, /* dimension 6 */ 4.724766, /* dimension 7 */ 4.058712, /* dimension 8 */ 3.298509, /* dimension 9 */ 2.550164, /* dimension 10 */ 1.884104, /* dimension 11 */ 1.335263, /* dimension 12 */ 0.910629, /* dimension 13 */ 0.599265, /* dimension 14 */ 0.381443, /* dimension 15 */ 0.235331, /* dimension 16 */ 0.140981, /* dimension 17 */ 0.082146, /* dimension 18 */ 0.046622, /* dimension 19 */ 0.025807, /* dimension 20 */ }; #if NUMDIMS > 20 # error "not enough precomputed sphere volumes" #endif #define UnitSphereVolume UnitSphereVolumes[NUMDIMS] #endif /*----------------------------------------------------------------------------- | Calculate the n-dimensional volume of the bounding sphere of a rectangle -----------------------------------------------------------------------------*/ #if 0 /* * A fast approximation to the volume of the bounding sphere for the * given Rect. By Paul B. */ RectReal RTreeRectSphericalVolume(struct Rect *R, struct RTree *t) { register struct Rect *r = R; register int i; RectReal maxsize = (RectReal) 0, c_size; assert(r); if (Undefined(r)) return (RectReal) 0; for (i = 0; i < t->ndims; i++) { c_size = r->boundary[i + NUMDIMS] - r->boundary[i]; if (c_size > maxsize) maxsize = c_size; } return (RectReal) (pow(maxsize / 2, NUMDIMS) * UnitSphereVolumes[t->ndims]); } #endif /* * The exact volume of the bounding sphere for the given Rect. */ RectReal RTreeRectSphericalVolume(struct Rect * r, struct RTree * t) { int i; double sum_of_squares = 0, radius, half_extent; assert(r); if (Undefined(r)) return (RectReal) 0; for (i = 0; i < t->ndims; i++) { half_extent = (r->boundary[i + NUMDIMS] - r->boundary[i]) / 2; sum_of_squares += half_extent * half_extent; } radius = sqrt(sum_of_squares); return (RectReal) (pow(radius, t->ndims) * UnitSphereVolumes[t->ndims]); } /*----------------------------------------------------------------------------- | Calculate the n-dimensional surface area of a rectangle -----------------------------------------------------------------------------*/ RectReal RTreeRectSurfaceArea(struct Rect * r, struct RTree * t) { int i, j; RectReal j_extent, sum = (RectReal) 0; assert(r); if (Undefined(r)) return (RectReal) 0; for (i = 0; i < t->ndims; i++) { RectReal face_area = (RectReal) 1; for (j = 0; j < t->ndims; j++) /* exclude i extent from product in this dimension */ if (i != j) { j_extent = r->boundary[j + NUMDIMS] - r->boundary[j]; face_area *= j_extent; } sum += face_area; } return 2 * sum; } /*----------------------------------------------------------------------------- | Calculate the n-dimensional margin of a rectangle | the margin is the sum of the lengths of the edges -----------------------------------------------------------------------------*/ RectReal RTreeRectMargin(struct Rect * r, struct RTree * t) { int i; RectReal margin = 0.0; assert(r); for (i = 0; i < t->ndims; i++) { margin += r->boundary[i + NUMDIMS] - r->boundary[i]; } return margin; } /*----------------------------------------------------------------------------- | Combine two rectangles, make one that includes both. -----------------------------------------------------------------------------*/ struct Rect RTreeCombineRect(struct Rect *r, struct Rect *rr, struct RTree *t) { int i, j; struct Rect new_rect; assert(r && rr); if (Undefined(r)) return *rr; if (Undefined(rr)) return *r; for (i = 0; i < t->ndims; i++) { new_rect.boundary[i] = MIN(r->boundary[i], rr->boundary[i]); j = i + NUMDIMS; new_rect.boundary[j] = MAX(r->boundary[j], rr->boundary[j]); } return new_rect; } /*----------------------------------------------------------------------------- | Decide whether two rectangles overlap. -----------------------------------------------------------------------------*/ int RTreeOverlap(struct Rect *r, struct Rect *s, struct RTree *t) { register int i, j; assert(r && s); for (i = 0; i < t->ndims; i++) { j = i + NUMDIMS; /* index for high sides */ if (r->boundary[i] > s->boundary[j] || s->boundary[i] > r->boundary[j]) { return FALSE; } } return TRUE; } /*----------------------------------------------------------------------------- | Decide whether rectangle r is contained in rectangle s. -----------------------------------------------------------------------------*/ int RTreeContained(struct Rect *r, struct Rect *s, struct RTree *t) { register int i, j, result; assert(r && s); /* same as in RTreeOverlap() */ /* undefined rect is contained in any other */ if (Undefined(r)) return TRUE; /* no rect (except an undefined one) is contained in an undef rect */ if (Undefined(s)) return FALSE; result = TRUE; for (i = 0; i < t->ndims; i++) { j = i + NUMDIMS; /* index for high sides */ result = result && r->boundary[i] >= s->boundary[i] && r->boundary[j] <= s->boundary[j]; } return result; }