GRASS Programmer's Manual
6.4.1(2011)
|
00001 00002 /**************************************************************************** 00003 * MODULE: R-Tree library 00004 * 00005 * AUTHOR(S): Antonin Guttman - original code 00006 * Daniel Green (green@superliminal.com) - major clean-up 00007 * and implementation of bounding spheres 00008 * 00009 * PURPOSE: Multidimensional index 00010 * 00011 * COPYRIGHT: (C) 2001 by the GRASS Development Team 00012 * 00013 * This program is free software under the GNU General Public 00014 * License (>=v2). Read the file COPYING that comes with GRASS 00015 * for details. 00016 *****************************************************************************/ 00017 00018 #include <stdio.h> 00019 #include <stdlib.h> 00020 #include "assert.h" 00021 #include "index.h" 00022 00023 #include <float.h> 00024 #include <math.h> 00025 #include <grass/gis.h> 00026 00027 #define BIG_NUM (FLT_MAX/4.0) 00028 00029 00030 #define Undefined(x) ((x)->boundary[0] > (x)->boundary[NUMDIMS]) 00031 #define MIN(a, b) ((a) < (b) ? (a) : (b)) 00032 #define MAX(a, b) ((a) > (b) ? (a) : (b)) 00033 00034 00035 /*----------------------------------------------------------------------------- 00036 | Initialize a rectangle to have all 0 coordinates. 00037 -----------------------------------------------------------------------------*/ 00038 void RTreeInitRect(struct Rect *R) 00039 { 00040 register struct Rect *r = R; 00041 register int i; 00042 00043 for (i = 0; i < NUMSIDES; i++) 00044 r->boundary[i] = (RectReal) 0; 00045 } 00046 00047 00048 /*----------------------------------------------------------------------------- 00049 | Return a rect whose first low side is higher than its opposite side - 00050 | interpreted as an undefined rect. 00051 -----------------------------------------------------------------------------*/ 00052 struct Rect RTreeNullRect(void) 00053 { 00054 struct Rect r; 00055 register int i; 00056 00057 r.boundary[0] = (RectReal) 1; 00058 r.boundary[NUMDIMS] = (RectReal) - 1; 00059 for (i = 1; i < NUMDIMS; i++) 00060 r.boundary[i] = r.boundary[i + NUMDIMS] = (RectReal) 0; 00061 return r; 00062 } 00063 00064 00065 #if 0 00066 00067 /*----------------------------------------------------------------------------- 00068 | Fills in random coordinates in a rectangle. 00069 | The low side is guaranteed to be less than the high side. 00070 -----------------------------------------------------------------------------*/ 00071 void RTreeRandomRect(struct Rect *R) 00072 { 00073 register struct Rect *r = R; 00074 register int i; 00075 register RectReal width; 00076 00077 for (i = 0; i < NUMDIMS; i++) { 00078 /* width from 1 to 1000 / 4, more small ones 00079 */ 00080 width = drand48() * (1000 / 4) + 1; 00081 00082 /* sprinkle a given size evenly but so they stay in [0,100] 00083 */ 00084 r->boundary[i] = drand48() * (1000 - width); /* low side */ 00085 r->boundary[i + NUMDIMS] = r->boundary[i] + width; /* high side */ 00086 } 00087 } 00088 00089 00090 /*----------------------------------------------------------------------------- 00091 | Fill in the boundaries for a random search rectangle. 00092 | Pass in a pointer to a rect that contains all the data, 00093 | and a pointer to the rect to be filled in. 00094 | Generated rect is centered randomly anywhere in the data area, 00095 | and has size from 0 to the size of the data area in each dimension, 00096 | i.e. search rect can stick out beyond data area. 00097 -----------------------------------------------------------------------------*/ 00098 void RTreeSearchRect(struct Rect *Search, struct Rect *Data) 00099 { 00100 register struct Rect *search = Search, *data = Data; 00101 register int i, j; 00102 register RectReal size, center; 00103 00104 assert(search); 00105 assert(data); 00106 00107 for (i = 0; i < NUMDIMS; i++) { 00108 j = i + NUMDIMS; /* index for high side boundary */ 00109 if (data->boundary[i] > -BIG_NUM && data->boundary[j] < BIG_NUM) { 00110 size = (drand48() * (data->boundary[j] - 00111 data->boundary[i] + 1)) / 2; 00112 center = data->boundary[i] + drand48() * 00113 (data->boundary[j] - data->boundary[i] + 1); 00114 search->boundary[i] = center - size / 2; 00115 search->boundary[j] = center + size / 2; 00116 } 00117 else { /* some open boundary, search entire dimension */ 00118 00119 search->boundary[i] = -BIG_NUM; 00120 search->boundary[j] = BIG_NUM; 00121 } 00122 } 00123 } 00124 00125 #endif 00126 00127 /*----------------------------------------------------------------------------- 00128 | Print out the data for a rectangle. 00129 -----------------------------------------------------------------------------*/ 00130 void RTreePrintRect(struct Rect *R, int depth) 00131 { 00132 register struct Rect *r = R; 00133 register int i; 00134 00135 assert(r); 00136 00137 RTreeTabIn(depth); 00138 fprintf(stdout, "rect:\n"); 00139 for (i = 0; i < NUMDIMS; i++) { 00140 RTreeTabIn(depth + 1); 00141 fprintf(stdout, "%f\t%f\n", r->boundary[i], r->boundary[i + NUMDIMS]); 00142 } 00143 } 00144 00145 /*----------------------------------------------------------------------------- 00146 | Calculate the n-dimensional volume of a rectangle 00147 -----------------------------------------------------------------------------*/ 00148 RectReal RTreeRectVolume(struct Rect *R) 00149 { 00150 register struct Rect *r = R; 00151 register int i; 00152 register RectReal volume = (RectReal) 1; 00153 00154 assert(r); 00155 if (Undefined(r)) 00156 return (RectReal) 0; 00157 00158 for (i = 0; i < NUMDIMS; i++) 00159 volume *= r->boundary[i + NUMDIMS] - r->boundary[i]; 00160 assert(volume >= 0.0); 00161 return volume; 00162 } 00163 00164 00165 /*----------------------------------------------------------------------------- 00166 | Define the NUMDIMS-dimensional volume the unit sphere in that dimension into 00167 | the symbol "UnitSphereVolume" 00168 | Note that if the gamma function is available in the math library and if the 00169 | compiler supports static initialization using functions, this is 00170 | easily computed for any dimension. If not, the value can be precomputed and 00171 | taken from a table. The following code can do it either way. 00172 -----------------------------------------------------------------------------*/ 00173 00174 #ifdef gamma 00175 00176 /* computes the volume of an N-dimensional sphere. */ 00177 /* derived from formule in "Regular Polytopes" by H.S.M Coxeter */ 00178 static double sphere_volume(double dimension) 00179 { 00180 double log_gamma, log_volume; 00181 00182 log_gamma = gamma(dimension / 2.0 + 1); 00183 log_volume = dimension / 2.0 * log(M_PI) - log_gamma; 00184 return exp(log_volume); 00185 } 00186 static const double UnitSphereVolume = sphere_volume(NUMDIMS); 00187 00188 #else 00189 00190 /* Precomputed volumes of the unit spheres for the first few dimensions */ 00191 const double UnitSphereVolumes[] = { 00192 0.000000, /* dimension 0 */ 00193 2.000000, /* dimension 1 */ 00194 3.141593, /* dimension 2 */ 00195 4.188790, /* dimension 3 */ 00196 4.934802, /* dimension 4 */ 00197 5.263789, /* dimension 5 */ 00198 5.167713, /* dimension 6 */ 00199 4.724766, /* dimension 7 */ 00200 4.058712, /* dimension 8 */ 00201 3.298509, /* dimension 9 */ 00202 2.550164, /* dimension 10 */ 00203 1.884104, /* dimension 11 */ 00204 1.335263, /* dimension 12 */ 00205 0.910629, /* dimension 13 */ 00206 0.599265, /* dimension 14 */ 00207 0.381443, /* dimension 15 */ 00208 0.235331, /* dimension 16 */ 00209 0.140981, /* dimension 17 */ 00210 0.082146, /* dimension 18 */ 00211 0.046622, /* dimension 19 */ 00212 0.025807, /* dimension 20 */ 00213 }; 00214 00215 #if NUMDIMS > 20 00216 # error "not enough precomputed sphere volumes" 00217 #endif 00218 #define UnitSphereVolume UnitSphereVolumes[NUMDIMS] 00219 00220 #endif 00221 00222 00223 /*----------------------------------------------------------------------------- 00224 | Calculate the n-dimensional volume of the bounding sphere of a rectangle 00225 -----------------------------------------------------------------------------*/ 00226 00227 #if 0 00228 /* 00229 * A fast approximation to the volume of the bounding sphere for the 00230 * given Rect. By Paul B. 00231 */ 00232 RectReal RTreeRectSphericalVolume(struct Rect *R) 00233 { 00234 register struct Rect *r = R; 00235 register int i; 00236 RectReal maxsize = (RectReal) 0, c_size; 00237 00238 assert(r); 00239 if (Undefined(r)) 00240 return (RectReal) 0; 00241 for (i = 0; i < NUMDIMS; i++) { 00242 c_size = r->boundary[i + NUMDIMS] - r->boundary[i]; 00243 if (c_size > maxsize) 00244 maxsize = c_size; 00245 } 00246 return (RectReal) (pow(maxsize / 2, NUMDIMS) * UnitSphereVolume); 00247 } 00248 #endif 00249 00250 /* 00251 * The exact volume of the bounding sphere for the given Rect. 00252 */ 00253 RectReal RTreeRectSphericalVolume(struct Rect * R) 00254 { 00255 register struct Rect *r = R; 00256 register int i; 00257 register double sum_of_squares = 0, radius; 00258 00259 assert(r); 00260 if (Undefined(r)) 00261 return (RectReal) 0; 00262 for (i = 0; i < NUMDIMS; i++) { 00263 double half_extent = (r->boundary[i + NUMDIMS] - r->boundary[i]) / 2; 00264 00265 sum_of_squares += half_extent * half_extent; 00266 } 00267 radius = sqrt(sum_of_squares); 00268 return (RectReal) (pow(radius, NUMDIMS) * UnitSphereVolume); 00269 } 00270 00271 00272 /*----------------------------------------------------------------------------- 00273 | Calculate the n-dimensional surface area of a rectangle 00274 -----------------------------------------------------------------------------*/ 00275 RectReal RTreeRectSurfaceArea(struct Rect * R) 00276 { 00277 register struct Rect *r = R; 00278 register int i, j; 00279 register RectReal sum = (RectReal) 0; 00280 00281 assert(r); 00282 if (Undefined(r)) 00283 return (RectReal) 0; 00284 00285 for (i = 0; i < NUMDIMS; i++) { 00286 RectReal face_area = (RectReal) 1; 00287 00288 for (j = 0; j < NUMDIMS; j++) 00289 /* exclude i extent from product in this dimension */ 00290 if (i != j) { 00291 RectReal j_extent = r->boundary[j + NUMDIMS] - r->boundary[j]; 00292 00293 face_area *= j_extent; 00294 } 00295 sum += face_area; 00296 } 00297 return 2 * sum; 00298 } 00299 00300 00301 00302 /*----------------------------------------------------------------------------- 00303 | Combine two rectangles, make one that includes both. 00304 -----------------------------------------------------------------------------*/ 00305 struct Rect RTreeCombineRect(struct Rect *R, struct Rect *Rr) 00306 { 00307 register struct Rect *r = R, *rr = Rr; 00308 register int i, j; 00309 struct Rect new_rect; 00310 00311 assert(r && rr); 00312 00313 if (Undefined(r)) 00314 return *rr; 00315 00316 if (Undefined(rr)) 00317 return *r; 00318 00319 for (i = 0; i < NUMDIMS; i++) { 00320 new_rect.boundary[i] = MIN(r->boundary[i], rr->boundary[i]); 00321 j = i + NUMDIMS; 00322 new_rect.boundary[j] = MAX(r->boundary[j], rr->boundary[j]); 00323 } 00324 return new_rect; 00325 } 00326 00327 00328 /*----------------------------------------------------------------------------- 00329 | Decide whether two rectangles overlap. 00330 -----------------------------------------------------------------------------*/ 00331 int RTreeOverlap(struct Rect *R, struct Rect *S) 00332 { 00333 register struct Rect *r = R, *s = S; 00334 register int i, j; 00335 00336 assert(r && s); 00337 00338 for (i = 0; i < NUMDIMS; i++) { 00339 j = i + NUMDIMS; /* index for high sides */ 00340 if (r->boundary[i] > s->boundary[j] || 00341 s->boundary[i] > r->boundary[j]) { 00342 return FALSE; 00343 } 00344 } 00345 return TRUE; 00346 } 00347 00348 00349 /*----------------------------------------------------------------------------- 00350 | Decide whether rectangle r is contained in rectangle s. 00351 -----------------------------------------------------------------------------*/ 00352 int RTreeContained(struct Rect *R, struct Rect *S) 00353 { 00354 register struct Rect *r = R, *s = S; 00355 register int i, j, result; 00356 00357 assert(r && s); /* same as in RTreeOverlap() */ 00358 00359 /* undefined rect is contained in any other */ 00360 if (Undefined(r)) 00361 return TRUE; 00362 00363 /* no rect (except an undefined one) is contained in an undef rect */ 00364 if (Undefined(s)) 00365 return FALSE; 00366 00367 result = TRUE; 00368 for (i = 0; i < NUMDIMS; i++) { 00369 j = i + NUMDIMS; /* index for high sides */ 00370 result = result && r->boundary[i] >= s->boundary[i] 00371 && r->boundary[j] <= s->boundary[j]; 00372 } 00373 return result; 00374 }