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