00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
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
00048
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
00067
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
00077
00078 width = drand48() * (1000 / 4) + 1;
00079
00080
00081
00082 r->boundary[i] = drand48() * (1000-width);
00083 r->boundary[i + NUMDIMS] = r->boundary[i] + width;
00084 }
00085 }
00086
00087
00088
00089
00090
00091
00092
00093
00094
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;
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
00119 {
00120 search->boundary[i] = -BIG_NUM;
00121 search->boundary[j] = BIG_NUM;
00122 }
00123 }
00124 }
00125
00126 #endif
00127
00128
00129
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
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
00167
00168
00169
00170
00171
00172
00173
00174 #ifdef gamma
00175
00176
00177
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
00190 const double UnitSphereVolumes[] = {
00191 0.000000,
00192 2.000000,
00193 3.141593,
00194 4.188790,
00195 4.934802,
00196 5.263789,
00197 5.167713,
00198 4.724766,
00199 4.058712,
00200 3.298509,
00201 2.550164,
00202 1.884104,
00203 1.335263,
00204 0.910629,
00205 0.599265,
00206 0.381443,
00207 0.235331,
00208 0.140981,
00209 0.082146,
00210 0.046622,
00211 0.025807,
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
00223
00224
00225 #if 0
00226
00227
00228
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
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
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
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
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
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;
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
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
00357 if (Undefined(r))
00358 return TRUE;
00359
00360
00361 if (Undefined(s))
00362 return FALSE;
00363
00364 result = TRUE;
00365 for (i = 0; i < NUMDIMS; i++)
00366 {
00367 j = i + NUMDIMS;
00368 result = result
00369 && r->boundary[i] >= s->boundary[i]
00370 && r->boundary[j] <= s->boundary[j];
00371 }
00372 return result;
00373 }