granicus.if.org Git - postgresql/blob - contrib/cube/cube.c

   1 /******************************************************************************
   2   This file contains routines that can be bound to a Postgres backend and
   3   called by the backend in the process of processing queries.  The calling
   4   format for these routines is dictated by Postgres architecture.
   5 ******************************************************************************/
   6
   7 #include "postgres.h"
   8
   9 #include <math.h>
  10
  11 #include "access/gist.h"
  12 #include "access/rtree.h"
  13 #include "utils/elog.h"
  14 #include "utils/palloc.h"
  15 #include "utils/builtins.h"
  16
  17 #include "cubedata.h"
  18
  19 #define max(a,b)                ((a) >  (b) ? (a) : (b))
  20 #define min(a,b)                ((a) <= (b) ? (a) : (b))
  21 #define abs(a)                  ((a) <  (0) ? (-a) : (a))
  22
  23 extern void set_parse_buffer(char *str);
  24 extern int      cube_yyparse();
  25
  26 /*
  27 ** Input/Output routines
  28 */
  29 NDBOX      *cube_in(char *str);
  30 char       *cube_out(NDBOX * cube);
  31
  32
  33 /*
  34 ** GiST support methods
  35 */
  36 bool            g_cube_consistent(GISTENTRY *entry, NDBOX * query, StrategyNumber strategy);
  37 GISTENTRY  *g_cube_compress(GISTENTRY *entry);
  38 GISTENTRY  *g_cube_decompress(GISTENTRY *entry);
  39 float      *g_cube_penalty(GISTENTRY *origentry, GISTENTRY *newentry, float *result);
  40 GIST_SPLITVEC *g_cube_picksplit(bytea *entryvec, GIST_SPLITVEC *v);
  41 bool            g_cube_leaf_consistent(NDBOX * key, NDBOX * query, StrategyNumber strategy);
  42 bool            g_cube_internal_consistent(NDBOX * key, NDBOX * query, StrategyNumber strategy);
  43 NDBOX      *g_cube_union(bytea *entryvec, int *sizep);
  44 NDBOX      *g_cube_binary_union(NDBOX * r1, NDBOX * r2, int *sizep);
  45 bool       *g_cube_same(NDBOX * b1, NDBOX * b2, bool *result);
  46
  47 /*
  48 ** R-tree suport functions
  49 */
  50 bool            cube_same(NDBOX * a, NDBOX * b);
  51 bool            cube_different(NDBOX * a, NDBOX * b);
  52 bool            cube_contains(NDBOX * a, NDBOX * b);
  53 bool            cube_contained(NDBOX * a, NDBOX * b);
  54 bool            cube_overlap(NDBOX * a, NDBOX * b);
  55 NDBOX      *cube_union(NDBOX * a, NDBOX * b);
  56 NDBOX      *cube_inter(NDBOX * a, NDBOX * b);
  57 float      *cube_size(NDBOX * a);
  58 void            rt_cube_size(NDBOX * a, float *sz);
  59
  60 /*
  61 ** These make no sense for this type, but R-tree wants them
  62 */
  63 bool            cube_over_left(NDBOX * a, NDBOX * b);
  64 bool            cube_over_right(NDBOX * a, NDBOX * b);
  65 bool            cube_left(NDBOX * a, NDBOX * b);
  66 bool            cube_right(NDBOX * a, NDBOX * b);
  67
  68 /*
  69 ** miscellaneous
  70 */
  71 bool            cube_lt(NDBOX * a, NDBOX * b);
  72 bool            cube_gt(NDBOX * a, NDBOX * b);
  73 float      *cube_distance(NDBOX * a, NDBOX * b);
  74
  75 /*
  76 ** Auxiliary funxtions
  77 */
  78 static float distance_1D(float a1, float a2, float b1, float b2);
  79 static NDBOX *swap_corners(NDBOX * a);
  80
  81
  82 /*****************************************************************************
  83  * Input/Output functions
  84  *****************************************************************************/
  85
  86 /* NdBox = [(lowerleft),(upperright)] */
  87 /* [(xLL(1)...xLL(N)),(xUR(1)...xUR(n))] */
  88 NDBOX *
  89 cube_in(char *str)
  90 {
  91         void       *result;
  92
  93         set_parse_buffer(str);
  94
  95         if (cube_yyparse(&result) != 0)
  96                 return NULL;
  97
  98         return ((NDBOX *) result);
  99 }
 100
 101 /*
 102  * You might have noticed a slight inconsistency between the following
 103  * declaration and the SQL definition:
 104  *         CREATE FUNCTION cube_out(opaque) RETURNS opaque ...
 105  * The reason is that the argument pass into cube_out is really just a
 106  * pointer. POSTGRES thinks all output functions are:
 107  *         char *out_func(char *);
 108  */
 109 char *
 110 cube_out(NDBOX * cube)
 111 {
 112         char       *result;
 113         char       *p;
 114         int                     equal = 1;
 115         int                     dim = cube->dim;
 116         int                     i;
 117
 118         if (cube == NULL)
 119                 return (NULL);
 120
 121         p = result = (char *) palloc(100);
 122
 123         /*
 124          * while printing the first (LL) corner, check if it is equal to the
 125          * scond one
 126          */
 127         p += sprintf(p, "(");
 128         for (i = 0; i < dim; i++)
 129         {
 130                 p += sprintf(p, "%g", cube->x[i]);
 131                 p += sprintf(p, ", ");
 132                 if (cube->x[i] != cube->x[i + dim])
 133                         equal = 0;
 134         }
 135         p -= 2;                                         /* get rid of the last ", " */
 136         p += sprintf(p, ")");
 137
 138         if (!equal)
 139         {
 140                 p += sprintf(p, ",(");
 141                 for (i = dim; i < dim * 2; i++)
 142                 {
 143                         p += sprintf(p, "%g", cube->x[i]);
 144                         p += sprintf(p, ", ");
 145                 }
 146                 p -= 2;
 147                 p += sprintf(p, ")");
 148         }
 149
 150         return (result);
 151 }
 152
 153
 154 /*****************************************************************************
 155  *                                                 GiST functions
 156  *****************************************************************************/
 157
 158 /*
 159 ** The GiST Consistent method for boxes
 160 ** Should return false if for all data items x below entry,
 161 ** the predicate x op query == FALSE, where op is the oper
 162 ** corresponding to strategy in the pg_amop table.
 163 */
 164 bool
 165 g_cube_consistent(GISTENTRY *entry,
 166                                   NDBOX * query,
 167                                   StrategyNumber strategy)
 168 {
 169
 170         /*
 171          * * if entry is not leaf, use g_cube_internal_consistent, * else use
 172          * g_cube_leaf_consistent
 173          */
 174         if (GIST_LEAF(entry))
 175                 return (g_cube_leaf_consistent((NDBOX *) (entry->pred), query, strategy));
 176         else
 177                 return (g_cube_internal_consistent((NDBOX *) (entry->pred), query, strategy));
 178 }
 179
 180
 181 /*
 182 ** The GiST Union method for boxes
 183 ** returns the minimal bounding box that encloses all the entries in entryvec
 184 */
 185 NDBOX *
 186 g_cube_union(bytea *entryvec, int *sizep)
 187 {
 188         int                     numranges,
 189                                 i;
 190         NDBOX      *out = (NDBOX *) NULL;
 191         NDBOX      *tmp;
 192
 193         /*
 194          * fprintf(stderr, "union\n");
 195          */
 196         numranges = (VARSIZE(entryvec) - VARHDRSZ) / sizeof(GISTENTRY);
 197         tmp = (NDBOX *) (((GISTENTRY *) (VARDATA(entryvec)))[0]).pred;
 198
 199         /*
 200          * sizep = sizeof(NDBOX); -- NDBOX has variable size
 201          */
 202         *sizep = tmp->size;
 203
 204         for (i = 1; i < numranges; i++)
 205         {
 206                 out = g_cube_binary_union(tmp, (NDBOX *)
 207                                                    (((GISTENTRY *) (VARDATA(entryvec)))[i]).pred,
 208                                                                   sizep);
 209
 210                 /*
 211                  * fprintf(stderr, "\t%s ^ %s -> %s\n", cube_out(tmp),
 212                  * cube_out((NDBOX *)(((GISTENTRY
 213                  * *)(VARDATA(entryvec)))[i]).pred), cube_out(out));
 214                  */
 215                 if (i > 1)
 216                         pfree(tmp);
 217                 tmp = out;
 218         }
 219
 220         return (out);
 221 }
 222
 223 /*
 224 ** GiST Compress and Decompress methods for boxes
 225 ** do not do anything.
 226 */
 227 GISTENTRY  *
 228 g_cube_compress(GISTENTRY *entry)
 229 {
 230         return (entry);
 231 }
 232
 233 GISTENTRY  *
 234 g_cube_decompress(GISTENTRY *entry)
 235 {
 236         return (entry);
 237 }
 238
 239 /*
 240 ** The GiST Penalty method for boxes
 241 ** As in the R-tree paper, we use change in area as our penalty metric
 242 */
 243 float *
 244 g_cube_penalty(GISTENTRY *origentry, GISTENTRY *newentry, float *result)
 245 {
 246         Datum           ud;
 247         float           tmp1,
 248                                 tmp2;
 249
 250         ud = (Datum) cube_union((NDBOX *) (origentry->pred), (NDBOX *) (newentry->pred));
 251         rt_cube_size((NDBOX *) ud, &tmp1);
 252         rt_cube_size((NDBOX *) (origentry->pred), &tmp2);
 253         *result = tmp1 - tmp2;
 254         pfree((char *) ud);
 255
 256         /*
 257          * fprintf(stderr, "penalty\n"); fprintf(stderr, "\t%g\n", *result);
 258          */
 259         return (result);
 260 }
 261
 262
 263
 264 /*
 265 ** The GiST PickSplit method for boxes
 266 ** We use Guttman's poly time split algorithm
 267 */
 268 GIST_SPLITVEC *
 269 g_cube_picksplit(bytea *entryvec,
 270                                  GIST_SPLITVEC *v)
 271 {
 272         OffsetNumber i,
 273                                 j;
 274         NDBOX      *datum_alpha,
 275                            *datum_beta;
 276         NDBOX      *datum_l,
 277                            *datum_r;
 278         NDBOX      *union_d,
 279                            *union_dl,
 280                            *union_dr;
 281         NDBOX      *inter_d;
 282         bool            firsttime;
 283         float           size_alpha,
 284                                 size_beta,
 285                                 size_union,
 286                                 size_inter;
 287         float           size_waste,
 288                                 waste;
 289         float           size_l,
 290                                 size_r;
 291         int                     nbytes;
 292         OffsetNumber seed_1 = 0,
 293                                 seed_2 = 0;
 294         OffsetNumber *left,
 295                            *right;
 296         OffsetNumber maxoff;
 297
 298         /*
 299          * fprintf(stderr, "picksplit\n");
 300          */
 301         maxoff = ((VARSIZE(entryvec) - VARHDRSZ) / sizeof(GISTENTRY)) - 2;
 302         nbytes = (maxoff + 2) * sizeof(OffsetNumber);
 303         v->spl_left = (OffsetNumber *) palloc(nbytes);
 304         v->spl_right = (OffsetNumber *) palloc(nbytes);
 305
 306         firsttime = true;
 307         waste = 0.0;
 308
 309         for (i = FirstOffsetNumber; i < maxoff; i = OffsetNumberNext(i))
 310         {
 311                 datum_alpha = (NDBOX *) (((GISTENTRY *) (VARDATA(entryvec)))[i].pred);
 312                 for (j = OffsetNumberNext(i); j <= maxoff; j = OffsetNumberNext(j))
 313                 {
 314                         datum_beta = (NDBOX *) (((GISTENTRY *) (VARDATA(entryvec)))[j].pred);
 315
 316                         /* compute the wasted space by unioning these guys */
 317                         /* size_waste = size_union - size_inter; */
 318                         union_d = (NDBOX *) cube_union(datum_alpha, datum_beta);
 319                         rt_cube_size(union_d, &size_union);
 320                         inter_d = (NDBOX *) cube_inter(datum_alpha, datum_beta);
 321                         rt_cube_size(inter_d, &size_inter);
 322                         size_waste = size_union - size_inter;
 323
 324                         pfree(union_d);
 325
 326                         if (inter_d != (NDBOX *) NULL)
 327                                 pfree(inter_d);
 328
 329                         /*
 330                          * are these a more promising split than what we've already
 331                          * seen?
 332                          */
 333
 334                         if (size_waste > waste || firsttime)
 335                         {
 336                                 waste = size_waste;
 337                                 seed_1 = i;
 338                                 seed_2 = j;
 339                                 firsttime = false;
 340                         }
 341                 }
 342         }
 343
 344         left = v->spl_left;
 345         v->spl_nleft = 0;
 346         right = v->spl_right;
 347         v->spl_nright = 0;
 348
 349         datum_alpha = (NDBOX *) (((GISTENTRY *) (VARDATA(entryvec)))[seed_1].pred);
 350         datum_l = (NDBOX *) cube_union(datum_alpha, datum_alpha);
 351         rt_cube_size((NDBOX *) datum_l, &size_l);
 352         datum_beta = (NDBOX *) (((GISTENTRY *) (VARDATA(entryvec)))[seed_2].pred);;
 353         datum_r = (NDBOX *) cube_union(datum_beta, datum_beta);
 354         rt_cube_size((NDBOX *) datum_r, &size_r);
 355
 356         /*
 357          * Now split up the regions between the two seeds.      An important
 358          * property of this split algorithm is that the split vector v has the
 359          * indices of items to be split in order in its left and right
 360          * vectors.  We exploit this property by doing a merge in the code
 361          * that actually splits the page.
 362          *
 363          * For efficiency, we also place the new index tuple in this loop. This
 364          * is handled at the very end, when we have placed all the existing
 365          * tuples and i == maxoff + 1.
 366          */
 367
 368         maxoff = OffsetNumberNext(maxoff);
 369         for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
 370         {
 371
 372                 /*
 373                  * If we've already decided where to place this item, just put it
 374                  * on the right list.  Otherwise, we need to figure out which page
 375                  * needs the least enlargement in order to store the item.
 376                  */
 377
 378                 if (i == seed_1)
 379                 {
 380                         *left++ = i;
 381                         v->spl_nleft++;
 382                         continue;
 383                 }
 384                 else if (i == seed_2)
 385                 {
 386                         *right++ = i;
 387                         v->spl_nright++;
 388                         continue;
 389                 }
 390
 391                 /* okay, which page needs least enlargement? */
 392                 datum_alpha = (NDBOX *) (((GISTENTRY *) (VARDATA(entryvec)))[i].pred);
 393                 union_dl = (NDBOX *) cube_union(datum_l, datum_alpha);
 394                 union_dr = (NDBOX *) cube_union(datum_r, datum_alpha);
 395                 rt_cube_size((NDBOX *) union_dl, &size_alpha);
 396                 rt_cube_size((NDBOX *) union_dr, &size_beta);
 397
 398                 /* pick which page to add it to */
 399                 if (size_alpha - size_l < size_beta - size_r)
 400                 {
 401                         pfree(datum_l);
 402                         pfree(union_dr);
 403                         datum_l = union_dl;
 404                         size_l = size_alpha;
 405                         *left++ = i;
 406                         v->spl_nleft++;
 407                 }
 408                 else
 409                 {
 410                         pfree(datum_r);
 411                         pfree(union_dl);
 412                         datum_r = union_dr;
 413                         size_r = size_alpha;
 414                         *right++ = i;
 415                         v->spl_nright++;
 416                 }
 417         }
 418         *left = *right = FirstOffsetNumber; /* sentinel value, see dosplit() */
 419
 420         v->spl_ldatum = (char *) datum_l;
 421         v->spl_rdatum = (char *) datum_r;
 422
 423         return v;
 424 }
 425
 426 /*
 427 ** Equality method
 428 */
 429 bool *
 430 g_cube_same(NDBOX * b1, NDBOX * b2, bool *result)
 431 {
 432         if (cube_same(b1, b2))
 433                 *result = TRUE;
 434         else
 435                 *result = FALSE;
 436
 437         /*
 438          * fprintf(stderr, "same: %s\n", (*result ? "TRUE" : "FALSE" ));
 439          */
 440         return (result);
 441 }
 442
 443 /*
 444 ** SUPPORT ROUTINES
 445 */
 446 bool
 447 g_cube_leaf_consistent(NDBOX * key,
 448                                            NDBOX * query,
 449                                            StrategyNumber strategy)
 450 {
 451         bool            retval;
 452
 453         /*
 454          * fprintf(stderr, "leaf_consistent, %d\n", strategy);
 455          */
 456         switch (strategy)
 457         {
 458                 case RTLeftStrategyNumber:
 459                         retval = (bool) cube_left(key, query);
 460                         break;
 461                 case RTOverLeftStrategyNumber:
 462                         retval = (bool) cube_over_left(key, query);
 463                         break;
 464                 case RTOverlapStrategyNumber:
 465                         retval = (bool) cube_overlap(key, query);
 466                         break;
 467                 case RTOverRightStrategyNumber:
 468                         retval = (bool) cube_over_right(key, query);
 469                         break;
 470                 case RTRightStrategyNumber:
 471                         retval = (bool) cube_right(key, query);
 472                         break;
 473                 case RTSameStrategyNumber:
 474                         retval = (bool) cube_same(key, query);
 475                         break;
 476                 case RTContainsStrategyNumber:
 477                         retval = (bool) cube_contains(key, query);
 478                         break;
 479                 case RTContainedByStrategyNumber:
 480                         retval = (bool) cube_contained(key, query);
 481                         break;
 482                 default:
 483                         retval = FALSE;
 484         }
 485         return (retval);
 486 }
 487
 488 bool
 489 g_cube_internal_consistent(NDBOX * key,
 490                                                    NDBOX * query,
 491                                                    StrategyNumber strategy)
 492 {
 493         bool            retval;
 494
 495         /*
 496          * fprintf(stderr, "internal_consistent, %d\n", strategy);
 497          */
 498         switch (strategy)
 499         {
 500                 case RTLeftStrategyNumber:
 501                 case RTOverLeftStrategyNumber:
 502                         retval = (bool) cube_over_left(key, query);
 503                         break;
 504                 case RTOverlapStrategyNumber:
 505                         retval = (bool) cube_overlap(key, query);
 506                         break;
 507                 case RTOverRightStrategyNumber:
 508                 case RTRightStrategyNumber:
 509                         retval = (bool) cube_right(key, query);
 510                         break;
 511                 case RTSameStrategyNumber:
 512                 case RTContainsStrategyNumber:
 513                         retval = (bool) cube_contains(key, query);
 514                         break;
 515                 case RTContainedByStrategyNumber:
 516                         retval = (bool) cube_overlap(key, query);
 517                         break;
 518                 default:
 519                         retval = FALSE;
 520         }
 521         return (retval);
 522 }
 523
 524 NDBOX *
 525 g_cube_binary_union(NDBOX * r1, NDBOX * r2, int *sizep)
 526 {
 527         NDBOX      *retval;
 528
 529         retval = cube_union(r1, r2);
 530         *sizep = retval->size;
 531
 532         return (retval);
 533 }
 534
 535
 536 /* cube_union */
 537 NDBOX *
 538 cube_union(NDBOX * box_a, NDBOX * box_b)
 539 {
 540         int                     i;
 541         NDBOX      *result;
 542         NDBOX      *a = swap_corners(box_a);
 543         NDBOX      *b = swap_corners(box_b);
 544
 545         if (a->dim >= b->dim)
 546         {
 547                 result = palloc(a->size);
 548                 result->size = a->size;
 549                 result->dim = a->dim;
 550         }
 551         else
 552         {
 553                 result = palloc(b->size);
 554                 result->size = b->size;
 555                 result->dim = b->dim;
 556         }
 557
 558         /* swap the box pointers if needed */
 559         if (a->dim < b->dim)
 560         {
 561                 NDBOX      *tmp = b;
 562
 563                 b = a;
 564                 a = tmp;
 565         }
 566
 567         /*
 568          * use the potentially smaller of the two boxes (b) to fill in the
 569          * result, padding absent dimensions with zeroes
 570          */
 571         for (i = 0; i < b->dim; i++)
 572         {
 573                 result->x[i] = b->x[i];
 574                 result->x[i + a->dim] = b->x[i + b->dim];
 575         }
 576         for (i = b->dim; i < a->dim; i++)
 577         {
 578                 result->x[i] = 0;
 579                 result->x[i + a->dim] = 0;
 580         }
 581
 582         /* compute the union */
 583         for (i = 0; i < a->dim; i++)
 584                 result->x[i] = min(a->x[i], result->x[i]);
 585         for (i = a->dim; i < a->dim * 2; i++)
 586                 result->x[i] = max(a->x[i], result->x[i]);
 587
 588         pfree(a);
 589         pfree(b);
 590
 591         return (result);
 592 }
 593
 594 /* cube_inter */
 595 NDBOX *
 596 cube_inter(NDBOX * box_a, NDBOX * box_b)
 597 {
 598         int                     i;
 599         NDBOX      *result;
 600         NDBOX      *a = swap_corners(box_a);
 601         NDBOX      *b = swap_corners(box_b);
 602
 603         if (a->dim >= b->dim)
 604         {
 605                 result = palloc(a->size);
 606                 result->size = a->size;
 607                 result->dim = a->dim;
 608         }
 609         else
 610         {
 611                 result = palloc(b->size);
 612                 result->size = b->size;
 613                 result->dim = b->dim;
 614         }
 615
 616         /* swap the box pointers if needed */
 617         if (a->dim < b->dim)
 618         {
 619                 NDBOX      *tmp = b;
 620
 621                 b = a;
 622                 a = tmp;
 623         }
 624
 625         /*
 626          * use the potentially  smaller of the two boxes (b) to fill in the
 627          * result, padding absent dimensions with zeroes
 628          */
 629         for (i = 0; i < b->dim; i++)
 630         {
 631                 result->x[i] = b->x[i];
 632                 result->x[i + a->dim] = b->x[i + b->dim];
 633         }
 634         for (i = b->dim; i < a->dim; i++)
 635         {
 636                 result->x[i] = 0;
 637                 result->x[i + a->dim] = 0;
 638         }
 639
 640         /* compute the intersection */
 641         for (i = 0; i < a->dim; i++)
 642                 result->x[i] = max(a->x[i], result->x[i]);
 643         for (i = a->dim; i < a->dim * 2; i++)
 644                 result->x[i] = min(a->x[i], result->x[i]);
 645
 646         pfree(a);
 647         pfree(b);
 648
 649         /*
 650          * Is it OK to return a non-null intersection for non-overlapping
 651          * boxes?
 652          */
 653         return (result);
 654 }
 655
 656 /* cube_size */
 657 float *
 658 cube_size(NDBOX * a)
 659 {
 660         int                     i,
 661                                 j;
 662         float      *result;
 663
 664         result = (float *) palloc(sizeof(float));
 665
 666         *result = 1.0;
 667         for (i = 0, j = a->dim; i < a->dim; i++, j++)
 668                 *result = (*result) * abs((a->x[j] - a->x[i]));
 669
 670         return (result);
 671 }
 672
 673 void
 674 rt_cube_size(NDBOX * a, float *size)
 675 {
 676         int                     i,
 677                                 j;
 678
 679         if (a == (NDBOX *) NULL)
 680                 *size = 0.0;
 681         else
 682         {
 683                 *size = 1.0;
 684                 for (i = 0, j = a->dim; i < a->dim; i++, j++)
 685                         *size = (*size) * abs((a->x[j] - a->x[i]));
 686         }
 687         return;
 688 }
 689
 690 /* The following four methods compare the projections of the boxes
 691    onto the 0-th coordinate axis. These methods are useless for dimensions
 692    larger than 2, but it seems that R-tree requires all its strategies
 693    map to real functions that return something */
 694
 695 /*      is the right edge of (a) located to the left of
 696         the right edge of (b)? */
 697 bool
 698 cube_over_left(NDBOX * box_a, NDBOX * box_b)
 699 {
 700         NDBOX      *a;
 701         NDBOX      *b;
 702
 703         if ((box_a == NULL) || (box_b == NULL))
 704                 return (FALSE);
 705
 706         a = swap_corners(box_a);
 707         b = swap_corners(box_b);
 708
 709         return (a->x[a->dim - 1] <= b->x[b->dim - 1] && !cube_left(a, b) && !cube_right(a, b));
 710 }
 711
 712 /*      is the left edge of (a) located to the right of
 713         the left edge of (b)? */
 714 bool
 715 cube_over_right(NDBOX * box_a, NDBOX * box_b)
 716 {
 717         NDBOX      *a;
 718         NDBOX      *b;
 719
 720         if ((box_a == NULL) || (box_b == NULL))
 721                 return (FALSE);
 722
 723         a = swap_corners(box_a);
 724         b = swap_corners(box_b);
 725
 726         return (a->x[a->dim - 1] >= b->x[b->dim - 1] && !cube_left(a, b) && !cube_right(a, b));
 727 }
 728
 729
 730 /* return 'true' if the projection of 'a' is
 731    entirely on the left of the projection of 'b' */
 732 bool
 733 cube_left(NDBOX * box_a, NDBOX * box_b)
 734 {
 735         NDBOX      *a;
 736         NDBOX      *b;
 737
 738         if ((box_a == NULL) || (box_b == NULL))
 739                 return (FALSE);
 740
 741         a = swap_corners(box_a);
 742         b = swap_corners(box_b);
 743
 744         return (a->x[a->dim - 1] < b->x[0]);
 745 }
 746
 747 /* return 'true' if the projection of 'a' is
 748    entirely on the right  of the projection of 'b' */
 749 bool
 750 cube_right(NDBOX * box_a, NDBOX * box_b)
 751 {
 752         NDBOX      *a;
 753         NDBOX      *b;
 754
 755         if ((box_a == NULL) || (box_b == NULL))
 756                 return (FALSE);
 757
 758         a = swap_corners(box_a);
 759         b = swap_corners(box_b);
 760
 761         return (a->x[0] > b->x[b->dim - 1]);
 762 }
 763
 764 /* make up a metric in which one box will be 'lower' than the other
 765    -- this can be useful for srting and to determine uniqueness */
 766 bool
 767 cube_lt(NDBOX * box_a, NDBOX * box_b)
 768 {
 769         int                     i;
 770         int                     dim;
 771         NDBOX      *a;
 772         NDBOX      *b;
 773
 774         if ((box_a == NULL) || (box_b == NULL))
 775                 return (FALSE);
 776
 777         a = swap_corners(box_a);
 778         b = swap_corners(box_b);
 779         dim = min(a->dim, b->dim);
 780
 781         /*
 782          * if all common dimensions are equal, the cube with more dimensions
 783          * wins
 784          */
 785         if (cube_same(a, b))
 786         {
 787                 if (a->dim < b->dim)
 788                         return (TRUE);
 789                 else
 790                         return (FALSE);
 791         }
 792
 793         /* compare the common dimensions */
 794         for (i = 0; i < dim; i++)
 795         {
 796                 if (a->x[i] > b->x[i])
 797                         return (FALSE);
 798                 if (a->x[i] < b->x[i])
 799                         return (TRUE);
 800         }
 801         for (i = 0; i < dim; i++)
 802         {
 803                 if (a->x[i + a->dim] > b->x[i + b->dim])
 804                         return (FALSE);
 805                 if (a->x[i + a->dim] < b->x[i + b->dim])
 806                         return (TRUE);
 807         }
 808
 809         /* compare extra dimensions to zero */
 810         if (a->dim > b->dim)
 811         {
 812                 for (i = dim; i < a->dim; i++)
 813                 {
 814                         if (a->x[i] > 0)
 815                                 return (FALSE);
 816                         if (a->x[i] < 0)
 817                                 return (TRUE);
 818                 }
 819                 for (i = 0; i < dim; i++)
 820                 {
 821                         if (a->x[i + a->dim] > 0)
 822                                 return (FALSE);
 823                         if (a->x[i + a->dim] < 0)
 824                                 return (TRUE);
 825                 }
 826         }
 827         if (a->dim < b->dim)
 828         {
 829                 for (i = dim; i < b->dim; i++)
 830                 {
 831                         if (b->x[i] > 0)
 832                                 return (TRUE);
 833                         if (b->x[i] < 0)
 834                                 return (FALSE);
 835                 }
 836                 for (i = 0; i < dim; i++)
 837                 {
 838                         if (b->x[i + b->dim] > 0)
 839                                 return (TRUE);
 840                         if (b->x[i + b->dim] < 0)
 841                                 return (FALSE);
 842                 }
 843         }
 844
 845         return (FALSE);
 846 }
 847
 848
 849 bool
 850 cube_gt(NDBOX * box_a, NDBOX * box_b)
 851 {
 852         int                     i;
 853         int                     dim;
 854         NDBOX      *a;
 855         NDBOX      *b;
 856
 857         if ((box_a == NULL) || (box_b == NULL))
 858                 return (FALSE);
 859
 860         a = swap_corners(box_a);
 861         b = swap_corners(box_b);
 862         dim = min(a->dim, b->dim);
 863
 864         /*
 865          * if all common dimensions are equal, the cube with more dimensions
 866          * wins
 867          */
 868         if (cube_same(a, b))
 869         {
 870                 if (a->dim > b->dim)
 871                         return (TRUE);
 872                 else
 873                         return (FALSE);
 874         }
 875
 876         /* compare the common dimensions */
 877         for (i = 0; i < dim; i++)
 878         {
 879                 if (a->x[i] < b->x[i])
 880                         return (FALSE);
 881                 if (a->x[i] > b->x[i])
 882                         return (TRUE);
 883         }
 884         for (i = 0; i < dim; i++)
 885         {
 886                 if (a->x[i + a->dim] < b->x[i + b->dim])
 887                         return (FALSE);
 888                 if (a->x[i + a->dim] > b->x[i + b->dim])
 889                         return (TRUE);
 890         }
 891
 892
 893         /* compare extra dimensions to zero */
 894         if (a->dim > b->dim)
 895         {
 896                 for (i = dim; i < a->dim; i++)
 897                 {
 898                         if (a->x[i] < 0)
 899                                 return (FALSE);
 900                         if (a->x[i] > 0)
 901                                 return (TRUE);
 902                 }
 903                 for (i = 0; i < dim; i++)
 904                 {
 905                         if (a->x[i + a->dim] < 0)
 906                                 return (FALSE);
 907                         if (a->x[i + a->dim] > 0)
 908                                 return (TRUE);
 909                 }
 910         }
 911         if (a->dim < b->dim)
 912         {
 913                 for (i = dim; i < b->dim; i++)
 914                 {
 915                         if (b->x[i] < 0)
 916                                 return (TRUE);
 917                         if (b->x[i] > 0)
 918                                 return (FALSE);
 919                 }
 920                 for (i = 0; i < dim; i++)
 921                 {
 922                         if (b->x[i + b->dim] < 0)
 923                                 return (TRUE);
 924                         if (b->x[i + b->dim] > 0)
 925                                 return (FALSE);
 926                 }
 927         }
 928
 929         return (FALSE);
 930 }
 931
 932
 933 /* Equal */
 934 bool
 935 cube_same(NDBOX * box_a, NDBOX * box_b)
 936 {
 937         int                     i;
 938         NDBOX      *a;
 939         NDBOX      *b;
 940
 941         if ((box_a == NULL) || (box_b == NULL))
 942                 return (FALSE);
 943
 944         a = swap_corners(box_a);
 945         b = swap_corners(box_b);
 946
 947         /* swap the box pointers if necessary */
 948         if (a->dim < b->dim)
 949         {
 950                 NDBOX      *tmp = b;
 951
 952                 b = a;
 953                 a = tmp;
 954         }
 955
 956         for (i = 0; i < b->dim; i++)
 957         {
 958                 if (a->x[i] != b->x[i])
 959                         return (FALSE);
 960                 if (a->x[i + a->dim] != b->x[i + b->dim])
 961                         return (FALSE);
 962         }
 963
 964         /*
 965          * all dimensions of (b) are compared to those of (a); instead of
 966          * those in (a) absent in (b), compare (a) to zero
 967          */
 968         for (i = b->dim; i < a->dim; i++)
 969         {
 970                 if (a->x[i] != 0)
 971                         return (FALSE);
 972                 if (a->x[i + a->dim] != 0)
 973                         return (FALSE);
 974         }
 975
 976         pfree(a);
 977         pfree(b);
 978
 979         return (TRUE);
 980 }
 981
 982 /* Different */
 983 bool
 984 cube_different(NDBOX * box_a, NDBOX * box_b)
 985 {
 986         return (!cube_same(box_a, box_b));
 987 }
 988
 989
 990 /* Contains */
 991 /* Box(A) CONTAINS Box(B) IFF pt(A) < pt(B) */
 992 bool
 993 cube_contains(NDBOX * box_a, NDBOX * box_b)
 994 {
 995         int                     i;
 996         NDBOX      *a;
 997         NDBOX      *b;
 998
 999         if ((box_a == NULL) || (box_b == NULL))
1000                 return (FALSE);
1001
1002         a = swap_corners(box_a);
1003         b = swap_corners(box_b);
1004
1005         if (a->dim < b->dim)
1006         {
1007
1008                 /*
1009                  * the further comparisons will make sense if the excess
1010                  * dimensions of (b) were zeroes
1011                  */
1012                 for (i = a->dim; i < b->dim; i++)
1013                 {
1014                         if (b->x[i] != 0)
1015                                 return (FALSE);
1016                         if (b->x[i + b->dim] != 0)
1017                                 return (FALSE);
1018                 }
1019         }
1020
1021         /* Can't care less about the excess dimensions of (a), if any */
1022         for (i = 0; i < min(a->dim, b->dim); i++)
1023         {
1024                 if (a->x[i] > b->x[i])
1025                         return (FALSE);
1026                 if (a->x[i + a->dim] < b->x[i + b->dim])
1027                         return (FALSE);
1028         }
1029
1030         pfree(a);
1031         pfree(b);
1032
1033         return (TRUE);
1034 }
1035
1036 /* Contained */
1037 /* Box(A) Contained by Box(B) IFF Box(B) Contains Box(A) */
1038 bool
1039 cube_contained(NDBOX * a, NDBOX * b)
1040 {
1041         if (cube_contains(b, a) == TRUE)
1042                 return (TRUE);
1043         else
1044                 return (FALSE);
1045 }
1046
1047 /* Overlap */
1048 /* Box(A) Overlap Box(B) IFF (pt(a)LL < pt(B)UR) && (pt(b)LL < pt(a)UR) */
1049 bool
1050 cube_overlap(NDBOX * box_a, NDBOX * box_b)
1051 {
1052         int                     i;
1053         NDBOX      *a;
1054         NDBOX      *b;
1055
1056         /*
1057          * This *very bad* error was found in the source: if ( (a==NULL) ||
1058          * (b=NULL) ) return(FALSE);
1059          */
1060         if ((box_a == NULL) || (box_b == NULL))
1061                 return (FALSE);
1062
1063         a = swap_corners(box_a);
1064         b = swap_corners(box_b);
1065
1066         /* swap the box pointers if needed */
1067         if (a->dim < b->dim)
1068         {
1069                 NDBOX      *tmp = b;
1070
1071                 b = a;
1072                 a = tmp;
1073         }
1074
1075         /* compare within the dimensions of (b) */
1076         for (i = 0; i < b->dim; i++)
1077         {
1078                 if (a->x[i] > b->x[i + b->dim])
1079                         return (FALSE);
1080                 if (a->x[i + a->dim] < b->x[i])
1081                         return (FALSE);
1082         }
1083
1084         /* compare to zero those dimensions in (a) absent in (b) */
1085         for (i = b->dim; i < a->dim; i++)
1086         {
1087                 if (a->x[i] > 0)
1088                         return (FALSE);
1089                 if (a->x[i + a->dim] < 0)
1090                         return (FALSE);
1091         }
1092
1093         pfree(a);
1094         pfree(b);
1095
1096         return (TRUE);
1097 }
1098
1099
1100 /* Distance */
1101 /* The distance is computed as a per axis sum of the squared distances
1102    between 1D projections of the boxes onto Cartesian axes. Assuming zero
1103    distance between overlapping projections, this metric coincides with the
1104    "common sense" geometric distance */
1105 float *
1106 cube_distance(NDBOX * a, NDBOX * b)
1107 {
1108         int                     i;
1109         double          d,
1110                                 distance;
1111         float      *result;
1112
1113         result = (float *) palloc(sizeof(float));
1114
1115         /* swap the box pointers if needed */
1116         if (a->dim < b->dim)
1117         {
1118                 NDBOX      *tmp = b;
1119
1120                 b = a;
1121                 a = tmp;
1122         }
1123
1124         distance = 0.0;
1125         /* compute within the dimensions of (b) */
1126         for (i = 0; i < b->dim; i++)
1127         {
1128                 d = distance_1D(a->x[i], a->x[i + a->dim], b->x[i], b->x[i + b->dim]);
1129                 distance += d * d;
1130         }
1131
1132         /* compute distance to zero for those dimensions in (a) absent in (b) */
1133         for (i = b->dim; i < a->dim; i++)
1134         {
1135                 d = distance_1D(a->x[i], a->x[i + a->dim], 0.0, 0.0);
1136                 distance += d * d;
1137         }
1138
1139         *result = (float) sqrt(distance);
1140
1141         return (result);
1142 }
1143
1144 static float
1145 distance_1D(float a1, float a2, float b1, float b2)
1146 {
1147         /* interval (a) is entirely on the left of (b) */
1148         if ((a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2))
1149                 return (min(b1, b2) - max(a1, a2));
1150
1151         /* interval (a) is entirely on the right of (b) */
1152         if ((a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2))
1153                 return (min(a1, a2) - max(b1, b2));
1154
1155         /* the rest are all sorts of intersections */
1156         return (0.0);
1157 }
1158
1159 /* normalize the box's co-ordinates by placing min(xLL,xUR) to LL
1160    and max(xLL,xUR) to UR
1161 */
1162 static NDBOX *
1163 swap_corners(NDBOX * a)
1164 {
1165         int                     i,
1166                                 j;
1167         NDBOX      *result;
1168
1169         result = palloc(a->size);
1170         result->size = a->size;
1171         result->dim = a->dim;
1172
1173         for (i = 0, j = a->dim; i < a->dim; i++, j++)
1174         {
1175                 result->x[i] = min(a->x[i], a->x[j]);
1176                 result->x[j] = max(a->x[i], a->x[j]);
1177         }
1178
1179         return (result);
1180 }