GRASS Programmer's Manual  6.4.1(2011)
rect.c
Go to the documentation of this file.
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 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines