rect.c

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

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5