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