src/rtree.c

   1 /* $Id$ */
   2
   3 /*
   4  *                            COPYRIGHT
   5  *
   6  *  PCB, interactive printed circuit board design
   7  *  Copyright (C) 1994,1995,1996 Thomas Nau
   8  *  Copyright (C) 1998,1999,2000,2001,2002,2003,2004 harry eaton
   9  *
  10  *  This program is free software; you can redistribute it and/or modify
  11  *  it under the terms of the GNU General Public License as published by
  12  *  the Free Software Foundation; either version 2 of the License, or
  13  *  (at your option) any later version.
  14  *
  15  *  This program is distributed in the hope that it will be useful,
  16  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  *  GNU General Public License for more details.
  19  *
  20  *  You should have received a copy of the GNU General Public License
  21  *  along with this program; if not, write to the Free Software
  22  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  23  *
  24  *  Contact addresses for paper mail and Email:
  25  *  harry eaton, 6697 Buttonhole Ct, Columbia, MD 21044 USA
  26  *  haceaton@aplcomm.jhuapl.edu
  27  *
  28  */
  29
  30 /* this file, rtree.c, was written and is
  31  * Copyright (c) 2004, harry eaton
  32  */
  33
  34 /* implements r-tree structures.
  35  * these should be much faster for the auto-router
  36  * because the recursive search is much more efficient
  37  * and that's where the auto-router spends all its time.
  38  */
  39 #ifdef HAVE_CONFIG_H
  40 #include "config.h"
  41 #endif
  42
  43 #include "global.h"
  44
  45 #include <assert.h>
  46 #include <setjmp.h>
  47
  48 #include "mymem.h"
  49
  50 #include "rtree.h"
  51
  52 #ifdef HAVE_LIBDMALLOC
  53 #include <dmalloc.h>
  54 #endif
  55
  56 RCSID ("$Id$");
  57
  58
  59 #define SLOW_ASSERTS
  60 /* All rectangles are closed on the bottom left and open on the
  61  * top right. i.e. they contain one corner point, but not the other.
  62  * This requires that the corner points not be equal!
  63  */
  64
  65 /* the number of entries in each rtree node
  66  * 4 - 7 seem to be pretty good settings
  67  */
  68 #define M_SIZE 6
  69
  70 /* it seems that sorting the leaf order slows us down
  71  * but sometimes gives better routes
  72  */
  73 #undef SORT
  74 #define SORT_NONLEAF
  75
  76 #define DELETE_BY_POINTER
  77
  78 typedef struct
  79 {
  80   const BoxType *bptr;          /* pointer to the box */
  81   BoxType bounds;               /* copy of the box for locality of reference */
  82 } Rentry;
  83
  84 struct rtree_node
  85 {
  86   BoxType box;                  /* bounds rectangle of this node */
  87   struct rtree_node *parent;    /* parent of this node, NULL = root */
  88   struct
  89   {
  90     unsigned is_leaf:1;         /* this is a leaf node */
  91     unsigned manage:31;         /* true==should free 'rect.bptr' if node is destroyed */
  92   }
  93   flags;
  94   union
  95   {
  96     struct rtree_node *kids[M_SIZE + 1];        /* when not leaf */
  97     Rentry rects[M_SIZE + 1];   /* when leaf */
  98   } u;
  99 };
 100
 101 #ifndef NDEBUG
 102 #ifdef SLOW_ASSERTS
 103 static int
 104 __r_node_is_good (struct rtree_node *node)
 105 {
 106   int i, flag = 1;
 107   int kind = -1;
 108   bool last = false;
 109
 110   if (node == NULL)
 111     return 1;
 112   for (i = 0; i < M_SIZE; i++)
 113     {
 114       if (node->flags.is_leaf)
 115         {
 116           if (!node->u.rects[i].bptr)
 117             {
 118               last = true;
 119               continue;
 120             }
 121           /* check that once one entry is empty, all the rest are too */
 122           if (node->u.rects[i].bptr && last)
 123             assert (0);
 124           /* check that the box makes sense */
 125           if (node->box.X1 > node->box.X2)
 126             assert (0);
 127           if (node->box.Y1 > node->box.Y2)
 128             assert (0);
 129           /* check that bounds is the same as the pointer */
 130           if (node->u.rects[i].bounds.X1 != node->u.rects[i].bptr->X1)
 131             assert (0);
 132           if (node->u.rects[i].bounds.Y1 != node->u.rects[i].bptr->Y1)
 133             assert (0);
 134           if (node->u.rects[i].bounds.X2 != node->u.rects[i].bptr->X2)
 135             assert (0);
 136           if (node->u.rects[i].bounds.Y2 != node->u.rects[i].bptr->Y2)
 137             assert (0);
 138           /* check that entries are within node bounds */
 139           if (node->u.rects[i].bounds.X1 < node->box.X1)
 140             assert (0);
 141           if (node->u.rects[i].bounds.X2 > node->box.X2)
 142             assert (0);
 143           if (node->u.rects[i].bounds.Y1 < node->box.Y1)
 144             assert (0);
 145           if (node->u.rects[i].bounds.Y2 > node->box.Y2)
 146             assert (0);
 147         }
 148       else
 149         {
 150           if (!node->u.kids[i])
 151             {
 152               last = true;
 153               continue;
 154             }
 155           /* make sure all children are the same type */
 156           if (kind == -1)
 157             kind = node->u.kids[i]->flags.is_leaf;
 158           else if (kind != node->u.kids[i]->flags.is_leaf)
 159             assert (0);
 160           /* check that once one entry is empty, all the rest are too */
 161           if (node->u.kids[i] && last)
 162             assert (0);
 163           /* check that entries are within node bounds */
 164           if (node->u.kids[i]->box.X1 < node->box.X1)
 165             assert (0);
 166           if (node->u.kids[i]->box.X2 > node->box.X2)
 167             assert (0);
 168           if (node->u.kids[i]->box.Y1 < node->box.Y1)
 169             assert (0);
 170           if (node->u.kids[i]->box.Y2 > node->box.Y2)
 171             assert (0);
 172         }
 173       flag <<= 1;
 174     }
 175   /* check that we're completely in the parent's bounds */
 176   if (node->parent)
 177     {
 178       if (node->parent->box.X1 > node->box.X1)
 179         assert (0);
 180       if (node->parent->box.X2 < node->box.X2)
 181         assert (0);
 182       if (node->parent->box.Y1 > node->box.Y1)
 183         assert (0);
 184       if (node->parent->box.Y2 < node->box.Y2)
 185         assert (0);
 186     }
 187   /* make sure overflow is empty */
 188   if (!node->flags.is_leaf && node->u.kids[i])
 189     assert (0);
 190   if (node->flags.is_leaf && node->u.rects[i].bptr)
 191     assert (0);
 192   return 1;
 193 }
 194
 195 /* check the whole tree from this node down for consistency */
 196 static bool
 197 __r_tree_is_good (struct rtree_node *node)
 198 {
 199   int i;
 200
 201   if (!node)
 202     return 1;
 203   if (!__r_node_is_good (node))
 204     assert (0);
 205   if (node->flags.is_leaf)
 206     return 1;
 207   for (i = 0; i < M_SIZE; i++)
 208     {
 209       if (!__r_tree_is_good (node->u.kids[i]))
 210         return 0;
 211     }
 212   return 1;
 213 }
 214 #endif
 215 #endif
 216
 217 #ifndef NDEBUG
 218 /* print out the tree */
 219 void
 220 __r_dump_tree (struct rtree_node *node, int depth)
 221 {
 222   int i, j;
 223   static int count;
 224   static double area;
 225
 226   if (depth == 0)
 227     {
 228       area = 0;
 229       count = 0;
 230     }
 231   area +=
 232     (node->box.X2 - node->box.X1) * (double) (node->box.Y2 - node->box.Y1);
 233   count++;
 234   for (i = 0; i < depth; i++)
 235     printf ("  ");
 236   if (!node->flags.is_leaf)
 237     {
 238       printf ("p=0x%p node X(%d, %d) Y(%d, %d)\n", (void *) node,
 239               node->box.X1, node->box.X2, node->box.Y1, node->box.Y2);
 240     }
 241   else
 242     {
 243       printf ("p=0x%p leaf manage(%02x) X(%d, %d) Y(%d, %d)\n", (void *) node,
 244               node->flags.manage, node->box.X1, node->box.X2, node->box.Y1,
 245               node->box.Y2);
 246       for (j = 0; j < M_SIZE; j++)
 247         {
 248           if (!node->u.rects[j].bptr)
 249             break;
 250           area +=
 251             (node->u.rects[j].bounds.X2 -
 252              node->u.rects[j].bounds.X1) *
 253             (double) (node->u.rects[j].bounds.Y2 -
 254                       node->u.rects[j].bounds.Y1);
 255           count++;
 256           for (i = 0; i < depth + 1; i++)
 257             printf ("  ");
 258           printf ("entry 0x%p X(%d, %d) Y(%d, %d)\n",
 259                   (void *) (node->u.rects[j].bptr),
 260                   node->u.rects[j].bounds.X1, node->u.rects[j].bounds.X2,
 261                   node->u.rects[j].bounds.Y1, node->u.rects[j].bounds.Y2);
 262         }
 263       return;
 264     }
 265   for (i = 0; i < M_SIZE; i++)
 266     if (node->u.kids[i])
 267       __r_dump_tree (node->u.kids[i], depth + 1);
 268   if (depth == 0)
 269     printf ("average box area is %g\n", area / count);
 270 }
 271 #endif
 272
 273 /* Sort the children or entries of a node
 274  * according to the largest side.
 275  */
 276 #ifdef SORT
 277 static int
 278 cmp_box (const BoxType * a, const BoxType * b)
 279 {
 280   /* compare two box coordinates so that the __r_search
 281    * will fail at the earliest comparison possible.
 282    * It needs to see the biggest X1 first, then the
 283    * smallest X2, the biggest Y1 and smallest Y2
 284    */
 285   if (a->X1 < b->X1)
 286     return 0;
 287   if (a->X1 > b->X1)
 288     return 1;
 289   if (a->X2 > b->X2)
 290     return 0;
 291   if (a->X2 < b->X2)
 292     return 1;
 293   if (a->Y1 < b->Y1)
 294     return 0;
 295   if (a->Y1 > b->Y1)
 296     return 1;
 297   if (a->Y2 > b->Y2)
 298     return 0;
 299   return 1;
 300 }
 301
 302 static void
 303 sort_node (struct rtree_node *node)
 304 {
 305   if (node->flags.is_leaf)
 306     {
 307       register Rentry *r, *i, temp;
 308
 309       for (r = &node->u.rects[1]; r->bptr; r++)
 310         {
 311           temp = *r;
 312           i = r - 1;
 313           while (i >= &node->u.rects[0])
 314             {
 315               if (cmp_box (&i->bounds, &r->bounds))
 316                 break;
 317               *(i + 1) = *i;
 318               i--;
 319             }
 320           *(i + 1) = temp;
 321         }
 322     }
 323 #ifdef SORT_NONLEAF
 324   else
 325     {
 326       register struct rtree_node **r, **i, *temp;
 327
 328       for (r = &node->u.kids[1]; *r; r++)
 329         {
 330           temp = *r;
 331           i = r - 1;
 332           while (i >= &node->u.kids[0])
 333             {
 334               if (cmp_box (&(*i)->box, &(*r)->box))
 335                 break;
 336               *(i + 1) = *i;
 337               i--;
 338             }
 339           *(i + 1) = temp;
 340         }
 341     }
 342 #endif
 343 }
 344 #else
 345 #define sort_node(x)
 346 #endif
 347
 348 /* set the node bounds large enough to encompass all
 349  * of the children's rectangles
 350  */
 351 static void
 352 adjust_bounds (struct rtree_node *node)
 353 {
 354   int i;
 355
 356   assert (node);
 357   assert (node->u.kids[0]);
 358   if (node->flags.is_leaf)
 359     {
 360       node->box = node->u.rects[0].bounds;
 361       for (i = 1; i < M_SIZE + 1; i++)
 362         {
 363           if (!node->u.rects[i].bptr)
 364             return;
 365           MAKEMIN (node->box.X1, node->u.rects[i].bounds.X1);
 366           MAKEMAX (node->box.X2, node->u.rects[i].bounds.X2);
 367           MAKEMIN (node->box.Y1, node->u.rects[i].bounds.Y1);
 368           MAKEMAX (node->box.Y2, node->u.rects[i].bounds.Y2);
 369         }
 370     }
 371   else
 372     {
 373       node->box = node->u.kids[0]->box;
 374       for (i = 1; i < M_SIZE + 1; i++)
 375         {
 376           if (!node->u.kids[i])
 377             return;
 378           MAKEMIN (node->box.X1, node->u.kids[i]->box.X1);
 379           MAKEMAX (node->box.X2, node->u.kids[i]->box.X2);
 380           MAKEMIN (node->box.Y1, node->u.kids[i]->box.Y1);
 381           MAKEMAX (node->box.Y2, node->u.kids[i]->box.Y2);
 382         }
 383     }
 384 }
 385
 386 /* create an r-tree from an unsorted list of boxes.
 387  * the r-tree will keep pointers into
 388  * it, so don't free the box list until you've called r_destroy_tree.
 389  * if you set 'manage' to true, r_destroy_tree will free your boxlist.
 390  */
 391 rtree_t *
 392 r_create_tree (const BoxType * boxlist[], int N, int manage)
 393 {
 394   rtree_t *rtree;
 395   struct rtree_node *node;
 396   int i;
 397
 398   assert (N >= 0);
 399   rtree = (rtree_t *)calloc (1, sizeof (*rtree));
 400   /* start with a single empty leaf node */
 401   node = (struct rtree_node *)calloc (1, sizeof (*node));
 402   node->flags.is_leaf = 1;
 403   node->parent = NULL;
 404   rtree->root = node;
 405   /* simple, just insert all of the boxes! */
 406   for (i = 0; i < N; i++)
 407     {
 408       assert (boxlist[i]);
 409       r_insert_entry (rtree, boxlist[i], manage);
 410     }
 411 #ifdef SLOW_ASSERTS
 412   assert (__r_tree_is_good (rtree->root));
 413 #endif
 414   return rtree;
 415 }
 416
 417 static void
 418 __r_destroy_tree (struct rtree_node *node)
 419 {
 420   int i, flag = 1;
 421
 422   if (node->flags.is_leaf)
 423     for (i = 0; i < M_SIZE; i++)
 424       {
 425         if (!node->u.rects[i].bptr)
 426           break;
 427         if (node->flags.manage & flag)
 428           free ((void *) node->u.rects[i].bptr);
 429         flag = flag << 1;
 430       }
 431   else
 432     for (i = 0; i < M_SIZE; i++)
 433       {
 434         if (!node->u.kids[i])
 435           break;
 436         __r_destroy_tree (node->u.kids[i]);
 437       }
 438   free (node);
 439 }
 440
 441 /* free the memory associated with an rtree. */
 442 void
 443 r_destroy_tree (rtree_t ** rtree)
 444 {
 445
 446   __r_destroy_tree ((*rtree)->root);
 447   free (*rtree);
 448   *rtree = NULL;
 449 }
 450
 451 typedef struct
 452 {
 453   int (*check_it) (const BoxType * region, void *cl);
 454   int (*found_it) (const BoxType * box, void *cl);
 455   void *closure;
 456 } r_arg;
 457
 458 /* most of the auto-routing time is spent in this routine
 459  * so some careful thought has been given to maximizing the speed
 460  *
 461  */
 462 int
 463 __r_search (struct rtree_node *node, const BoxType * query, r_arg * arg)
 464 {
 465   assert (node);
 466   /** assert that starting_region is well formed */
 467   assert (query->X1 < query->X2 && query->Y1 < query->Y2);
 468   assert (node->box.X1 < query->X2 && node->box.X2 > query->X1 &&
 469           node->box.Y1 < query->Y2 && node->box.Y2 > query->Y1);
 470 #ifdef SLOW_ASSERTS
 471   /** assert that node is well formed */
 472   assert (__r_node_is_good (node));
 473 #endif
 474   /* the check for bounds is done before entry. This saves the overhead
 475    * of building/destroying the stack frame for each bounds that fails
 476    * to intersect, which is the most common condition.
 477    */
 478   if (node->flags.is_leaf)
 479     {
 480       register int i;
 481       if (arg->found_it)        /* test this once outside of loop */
 482         {
 483           register int seen = 0;
 484           for (i = 0; node->u.rects[i].bptr; i++)
 485             {
 486               if ((node->u.rects[i].bounds.X1 < query->X2) &&
 487                   (node->u.rects[i].bounds.X2 > query->X1) &&
 488                   (node->u.rects[i].bounds.Y1 < query->Y2) &&
 489                   (node->u.rects[i].bounds.Y2 > query->Y1) &&
 490                   arg->found_it (node->u.rects[i].bptr, arg->closure))
 491                 seen++;
 492             }
 493           return seen;
 494         }
 495       else
 496         {
 497           register int seen = 0;
 498           for (i = 0; node->u.rects[i].bptr; i++)
 499             {
 500               if ((node->u.rects[i].bounds.X1 < query->X2) &&
 501                   (node->u.rects[i].bounds.X2 > query->X1) &&
 502                   (node->u.rects[i].bounds.Y1 < query->Y2) &&
 503                   (node->u.rects[i].bounds.Y2 > query->Y1))
 504                 seen++;
 505             }
 506           return seen;
 507         }
 508     }
 509
 510   /* not a leaf, recurse on lower nodes */
 511   if (arg->check_it)
 512     {
 513       int seen = 0;
 514       struct rtree_node **n;
 515       for (n = &node->u.kids[0]; *n; n++)
 516         {
 517           if ((*n)->box.X1 >= query->X2 ||
 518               (*n)->box.X2 <= query->X1 ||
 519               (*n)->box.Y1 >= query->Y2 || (*n)->box.Y2 <= query->Y1 ||
 520               !arg->check_it (&(*n)->box, arg->closure))
 521             continue;
 522           seen += __r_search (*n, query, arg);
 523         }
 524       return seen;
 525     }
 526   else
 527     {
 528       int seen = 0;
 529       struct rtree_node **n;
 530       for (n = &node->u.kids[0]; *n; n++)
 531         {
 532           if ((*n)->box.X1 >= query->X2 ||
 533               (*n)->box.X2 <= query->X1 ||
 534               (*n)->box.Y1 >= query->Y2 || (*n)->box.Y2 <= query->Y1)
 535             continue;
 536           seen += __r_search (*n, query, arg);
 537         }
 538       return seen;
 539     }
 540 }
 541
 542 /* Parameterized search in the rtree.
 543  * Returns the number of rectangles found.
 544  * calls found_rectangle for each intersection seen
 545  * and calls check_region with the current sub-region
 546  * to see whether deeper searching is desired
 547  */
 548 int
 549 r_search (rtree_t * rtree, const BoxType * query,
 550           int (*check_region) (const BoxType * region, void *cl),
 551           int (*found_rectangle) (const BoxType * box, void *cl), void *cl)
 552 {
 553   r_arg arg;
 554
 555   if (!rtree || rtree->size < 1)
 556     return 0;
 557   if (query)
 558     {
 559 #ifdef SLOW_ASSERTS
 560   assert (__r_tree_is_good (rtree->root));
 561 #endif
 562 #ifdef DEBUG
 563       if (query->X2 <= query->X1 || query->Y2 <= query->Y1)
 564         return 0;
 565 #endif
 566       /* check this box. If it's not touched we're done here */
 567       if (rtree->root->box.X1 >= query->X2 ||
 568           rtree->root->box.X2 <= query->X1 ||
 569           rtree->root->box.Y1 >= query->Y2
 570           || rtree->root->box.Y2 <= query->Y1)
 571         return 0;
 572       arg.check_it = check_region;
 573       arg.found_it = found_rectangle;
 574       arg.closure = cl;
 575       return __r_search (rtree->root, query, &arg);
 576     }
 577   else
 578     {
 579       arg.check_it = check_region;
 580       arg.found_it = found_rectangle;
 581       arg.closure = cl;
 582       return __r_search (rtree->root, &rtree->root->box, &arg);
 583     }
 584 }
 585
 586 /*------ r_region_is_empty ------*/
 587 static int
 588 __r_region_is_empty_rect_in_reg (const BoxType * box, void *cl)
 589 {
 590   jmp_buf *envp = (jmp_buf *) cl;
 591   longjmp (*envp, 1);           /* found one! */
 592 }
 593
 594 /* return 0 if there are any rectangles in the given region. */
 595 int
 596 r_region_is_empty (rtree_t * rtree, const BoxType * region)
 597 {
 598   jmp_buf env;
 599 #ifndef NDEBUG
 600   int r;
 601 #endif
 602
 603   if (setjmp (env))
 604     return 0;
 605 #ifndef NDEBUG
 606   r =
 607 #endif
 608     r_search (rtree, region, NULL, __r_region_is_empty_rect_in_reg, &env);
 609 #ifndef NDEBUG
 610   assert (r == 0);              /* otherwise longjmp would have been called */
 611 #endif
 612   return 1;                     /* no rectangles found */
 613 }
 614
 615 struct centroid
 616 {
 617   float x, y, area;
 618 };
 619
 620 /* split the node into two nodes putting clusters in each
 621  * use the k-means clustering algorithm
 622  */
 623 struct rtree_node *
 624 find_clusters (struct rtree_node *node)
 625 {
 626   float total_a, total_b;
 627   float a_X, a_Y, b_X, b_Y;
 628   bool belong[M_SIZE + 1];
 629   struct centroid center[M_SIZE + 1];
 630   int clust_a, clust_b, tries;
 631   int a_manage = 0, b_manage = 0;
 632   int i, old_ax, old_ay, old_bx, old_by;
 633   struct rtree_node *new_node;
 634   BoxType *b;
 635
 636   for (i = 0; i < M_SIZE + 1; i++)
 637     {
 638       if (node->flags.is_leaf)
 639         b = &(node->u.rects[i].bounds);
 640       else
 641         b = &(node->u.kids[i]->box);
 642       center[i].x = 0.5 * (b->X1 + b->X2);
 643       center[i].y = 0.5 * (b->Y1 + b->Y2);
 644       /* adding 1 prevents zero area */
 645       center[i].area = 1. + (float) (b->X2 - b->X1) * (float) (b->Y2 - b->Y1);
 646     }
 647   /* starting 'A' cluster center */
 648   a_X = center[0].x;
 649   a_Y = center[0].y;
 650   /* starting 'B' cluster center */
 651   b_X = center[M_SIZE].x;
 652   b_Y = center[M_SIZE].y;
 653   /* don't allow the same cluster centers */
 654   if (b_X == a_X && b_Y == a_Y)
 655     {
 656       b_X += 10000;
 657       a_Y -= 10000;
 658     }
 659   for (tries = 0; tries < M_SIZE; tries++)
 660     {
 661       old_ax = (int) a_X;
 662       old_ay = (int) a_Y;
 663       old_bx = (int) b_X;
 664       old_by = (int) b_Y;
 665       clust_a = clust_b = 0;
 666       for (i = 0; i < M_SIZE + 1; i++)
 667         {
 668           float dist1, dist2;
 669
 670           dist1 = SQUARE (a_X - center[i].x) + SQUARE (a_Y - center[i].y);
 671           dist2 = SQUARE (b_X - center[i].x) + SQUARE (b_Y - center[i].y);
 672           if (dist1 * (clust_a + M_SIZE / 2) < dist2 * (clust_b + M_SIZE / 2))
 673             {
 674               belong[i] = true;
 675               clust_a++;
 676             }
 677           else
 678             {
 679               belong[i] = false;
 680               clust_b++;
 681             }
 682         }
 683       /* kludge to fix degenerate cases */
 684       if (clust_a == M_SIZE + 1)
 685         belong[M_SIZE / 2] = false;
 686       else if (clust_b == M_SIZE + 1)
 687         belong[M_SIZE / 2] = true;
 688       /* compute new center of gravity of clusters */
 689       total_a = total_b = 0;
 690       a_X = a_Y = b_X = b_Y = 0;
 691       for (i = 0; i < M_SIZE + 1; i++)
 692         {
 693           if (belong[i])
 694             {
 695               a_X += center[i].x * center[i].area;
 696               a_Y += center[i].y * center[i].area;
 697               total_a += center[i].area;
 698             }
 699           else
 700             {
 701               b_X += center[i].x * center[i].area;
 702               b_Y += center[i].y * center[i].area;
 703               total_b += center[i].area;
 704             }
 705         }
 706       a_X /= total_a;
 707       a_Y /= total_a;
 708       b_X /= total_b;
 709       b_Y /= total_b;
 710       if (old_ax == (int) a_X && old_ay == (int) a_Y &&
 711           old_bx == (int) b_X && old_by == (int) b_Y)
 712         break;
 713     }
 714   /* Now 'belong' has the partition map */
 715   new_node = (struct rtree_node *)calloc (1, sizeof (*new_node));
 716   new_node->parent = node->parent;
 717   new_node->flags.is_leaf = node->flags.is_leaf;
 718   clust_a = clust_b = 0;
 719   if (node->flags.is_leaf)
 720     {
 721       int flag, a_flag, b_flag;
 722       flag = a_flag = b_flag = 1;
 723       for (i = 0; i < M_SIZE + 1; i++)
 724         {
 725           if (belong[i])
 726             {
 727               node->u.rects[clust_a++] = node->u.rects[i];
 728               if (node->flags.manage & flag)
 729                 a_manage |= a_flag;
 730               a_flag <<= 1;
 731             }
 732           else
 733             {
 734               new_node->u.rects[clust_b++] = node->u.rects[i];
 735               if (node->flags.manage & flag)
 736                 b_manage |= b_flag;
 737               b_flag <<= 1;
 738             }
 739           flag <<= 1;
 740         }
 741     }
 742   else
 743     {
 744       for (i = 0; i < M_SIZE + 1; i++)
 745         {
 746           if (belong[i])
 747             node->u.kids[clust_a++] = node->u.kids[i];
 748           else
 749             {
 750               node->u.kids[i]->parent = new_node;
 751               new_node->u.kids[clust_b++] = node->u.kids[i];
 752             }
 753         }
 754     }
 755   node->flags.manage = a_manage;
 756   new_node->flags.manage = b_manage;
 757   assert (clust_a != 0);
 758   assert (clust_b != 0);
 759   if (node->flags.is_leaf)
 760     for (; clust_a < M_SIZE + 1; clust_a++)
 761       node->u.rects[clust_a].bptr = NULL;
 762   else
 763     for (; clust_a < M_SIZE + 1; clust_a++)
 764       node->u.kids[clust_a] = NULL;
 765   adjust_bounds (node);
 766   sort_node (node);
 767   adjust_bounds (new_node);
 768   sort_node (new_node);
 769   return (new_node);
 770 }
 771
 772 /* split a node according to clusters
 773  */
 774 static void
 775 split_node (struct rtree_node *node)
 776 {
 777   int i;
 778   struct rtree_node *new_node;
 779
 780   assert (node);
 781   assert (node->flags.is_leaf ? (void *) node->u.rects[M_SIZE].
 782           bptr : (void *) node->u.kids[M_SIZE]);
 783   new_node = find_clusters (node);
 784   if (node->parent == NULL)     /* split root node */
 785     {
 786       struct rtree_node *second;
 787
 788       second = (struct rtree_node *)calloc (1, sizeof (*second));
 789       *second = *node;
 790       if (!second->flags.is_leaf)
 791         for (i = 0; i < M_SIZE; i++)
 792           if (second->u.kids[i])
 793             second->u.kids[i]->parent = second;
 794       node->flags.is_leaf = 0;
 795       node->flags.manage = 0;
 796       second->parent = new_node->parent = node;
 797       node->u.kids[0] = new_node;
 798       node->u.kids[1] = second;
 799       for (i = 2; i < M_SIZE + 1; i++)
 800         node->u.kids[i] = NULL;
 801       adjust_bounds (node);
 802       sort_node (node);
 803 #ifdef SLOW_ASSERTS
 804       assert (__r_tree_is_good (node));
 805 #endif
 806       return;
 807     }
 808   for (i = 0; i < M_SIZE; i++)
 809     if (!node->parent->u.kids[i])
 810       break;
 811   node->parent->u.kids[i] = new_node;
 812 #ifdef SLOW_ASSERTS
 813   assert (__r_node_is_good (node));
 814   assert (__r_node_is_good (new_node));
 815 #endif
 816   if (i < M_SIZE)
 817     {
 818 #ifdef SLOW_ASSERTS
 819       assert (__r_node_is_good (node->parent));
 820 #endif
 821       sort_node (node->parent);
 822       return;
 823     }
 824   split_node (node->parent);
 825 }
 826
 827 static inline int
 828 contained (struct rtree_node *node, const BoxType * query)
 829 {
 830   if (node->box.X1 > query->X1 || node->box.X2 < query->X2 ||
 831       node->box.Y1 > query->Y1 || node->box.Y2 < query->Y2)
 832     return 0;
 833   return 1;
 834 }
 835
 836
 837 static inline double
 838 penalty (struct rtree_node *node, const BoxType * query)
 839 {
 840   double score;
 841
 842   /* Compute the area penalty for inserting here and return.
 843    * The penalty is the increase in area necessary
 844    * to include the query->
 845    */
 846   score  = (MAX (node->box.X2, query->X2) - MIN (node->box.X1, query->X1));
 847   score *= (MAX (node->box.Y2, query->Y2) - MIN (node->box.Y1, query->Y1));
 848   score -=
 849     (double)(node->box.X2 - node->box.X1) *
 850     (double)(node->box.Y2 - node->box.Y1);
 851   return score;
 852 }
 853
 854 static void
 855 __r_insert_node (struct rtree_node *node, const BoxType * query,
 856                  int manage, bool force)
 857 {
 858
 859 #ifdef SLOW_ASSERTS
 860   assert (__r_node_is_good (node));
 861 #endif
 862   /* Ok, at this point we must already enclose the query or we're forcing
 863    * this node to expand to enclose it, so if we're a leaf, simply store
 864    * the query here
 865    */
 866
 867   if (node->flags.is_leaf)
 868     {
 869       register int i;
 870
 871       if (UNLIKELY (manage))
 872         {
 873           register int flag = 1;
 874
 875           for (i = 0; i < M_SIZE; i++)
 876             {
 877               if (!node->u.rects[i].bptr)
 878                 break;
 879               flag <<= 1;
 880             }
 881           node->flags.manage |= flag;
 882         }
 883       else
 884         {
 885           for (i = 0; i < M_SIZE; i++)
 886             if (!node->u.rects[i].bptr)
 887               break;
 888         }
 889       /* the node always has an extra space available */
 890       node->u.rects[i].bptr = query;
 891       node->u.rects[i].bounds = *query;
 892       /* first entry in node determines initial bounding box */
 893       if (i == 0)
 894         node->box = *query;
 895       else if (force)
 896         {
 897           MAKEMIN (node->box.X1, query->X1);
 898           MAKEMAX (node->box.X2, query->X2);
 899           MAKEMIN (node->box.Y1, query->Y1);
 900           MAKEMAX (node->box.Y2, query->Y2);
 901         }
 902       if (i < M_SIZE)
 903         {
 904           sort_node (node);
 905           return;
 906         }
 907       /* we must split the node */
 908       split_node (node);
 909       return;
 910     }
 911   else
 912     {
 913       int i;
 914       struct rtree_node *best_node;
 915       double score, best_score;
 916
 917       if (force)
 918         {
 919           MAKEMIN (node->box.X1, query->X1);
 920           MAKEMAX (node->box.X2, query->X2);
 921           MAKEMIN (node->box.Y1, query->Y1);
 922           MAKEMAX (node->box.Y2, query->Y2);
 923         }
 924
 925       /* this node encloses it, but it's not a leaf, so descend the tree */
 926
 927       /* First check if any children actually encloses it */
 928       assert (node->u.kids[0]);
 929       for (i = 0; i < M_SIZE; i++)
 930         {
 931           if (!node->u.kids[i])
 932             break;
 933           if (contained (node->u.kids[i], query))
 934             {
 935               __r_insert_node (node->u.kids[i], query, manage, false);
 936               sort_node (node);
 937               return;
 938             }
 939         }
 940
 941       /* see if there is room for a new leaf node */
 942       if (node->u.kids[0]->flags.is_leaf && i < M_SIZE)
 943         {
 944           struct rtree_node *new_node;
 945           new_node = (struct rtree_node *)calloc (1, sizeof (*new_node));
 946           new_node->parent = node;
 947           new_node->flags.is_leaf = true;
 948           node->u.kids[i] = new_node;
 949           new_node->u.rects[0].bptr = query;
 950           new_node->u.rects[0].bounds = *query;
 951           new_node->box = *query;
 952           if (UNLIKELY (manage))
 953             new_node->flags.manage = 1;
 954           sort_node (node);
 955           return;
 956         }
 957
 958       /* Ok, so we're still here - look for the best child to push it into */
 959       best_score = penalty (node->u.kids[0], query);
 960       best_node = node->u.kids[0];
 961       for (i = 1; i < M_SIZE; i++)
 962         {
 963           if (!node->u.kids[i])
 964             break;
 965           score = penalty (node->u.kids[i], query);
 966           if (score < best_score)
 967             {
 968               best_score = score;
 969               best_node = node->u.kids[i];
 970             }
 971         }
 972       __r_insert_node (best_node, query, manage, true);
 973       sort_node (node);
 974       return;
 975     }
 976 }
 977
 978 void
 979 r_insert_entry (rtree_t * rtree, const BoxType * which, int man)
 980 {
 981   assert (which);
 982   assert (which->X1 <= which->X2);
 983   assert (which->Y1 <= which->Y2);
 984   /* recursively search the tree for the best leaf node */
 985   assert (rtree->root);
 986   __r_insert_node (rtree->root, which, man,
 987                    rtree->root->box.X1 > which->X1
 988                    || rtree->root->box.X2 < which->X2
 989                    || rtree->root->box.Y1 > which->Y1
 990                    || rtree->root->box.Y2 < which->Y2);
 991   rtree->size++;
 992 }
 993
 994 bool
 995 __r_delete (struct rtree_node *node, const BoxType * query)
 996 {
 997   int i, flag, mask, a;
 998
 999   /* the tree might be inconsistent during delete */
1000   if (query->X1 < node->box.X1 || query->Y1 < node->box.Y1
1001       || query->X2 > node->box.X2 || query->Y2 > node->box.Y2)
1002     return false;
1003   if (!node->flags.is_leaf)
1004     {
1005       for (i = 0; i < M_SIZE; i++)
1006         {
1007           /* if this is us being removed, free and copy over */
1008           if (node->u.kids[i] == (struct rtree_node *) query)
1009             {
1010               free ((void *) query);
1011               for (; i < M_SIZE; i++)
1012                 {
1013                   node->u.kids[i] = node->u.kids[i + 1];
1014                   if (!node->u.kids[i])
1015                     break;
1016                 }
1017               /* nobody home here now ? */
1018               if (!node->u.kids[0])
1019                 {
1020                   if (!node->parent)
1021                     {
1022                       /* wow, the root is empty! */
1023                       node->flags.is_leaf = 1;
1024                       /* changing type of node, be sure it's all zero */
1025                       for (i = 1; i < M_SIZE + 1; i++)
1026                         node->u.rects[i].bptr = NULL;
1027                       return true;
1028                     }
1029                   return (__r_delete (node->parent, &node->box));
1030                 }
1031               else
1032                 /* propegate boundary adjust upward */
1033                 while (node)
1034                   {
1035                     adjust_bounds (node);
1036                     node = node->parent;
1037                   }
1038               return true;
1039             }
1040           if (node->u.kids[i])
1041             {
1042               if (__r_delete (node->u.kids[i], query))
1043                 return true;
1044             }
1045           else
1046             break;
1047         }
1048       return false;
1049     }
1050   /* leaf node here */
1051   mask = 0;
1052   a = 1;
1053   for (i = 0; i < M_SIZE; i++)
1054     {
1055 #ifdef DELETE_BY_POINTER
1056       if (!node->u.rects[i].bptr || node->u.rects[i].bptr == query)
1057 #else
1058       if (node->u.rects[i].bounds.X1 == query->X1 &&
1059           node->u.rects[i].bounds.X2 == query->X2 &&
1060           node->u.rects[i].bounds.Y1 == query->Y1 &&
1061           node->u.rects[i].bounds.Y2 == query->Y2)
1062 #endif
1063         break;
1064       mask |= a;
1065       a <<= 1;
1066     }
1067   if (!node->u.rects[i].bptr)
1068     return false;               /* not at this leaf */
1069   if (node->flags.manage & a)
1070     {
1071       free ((void *) node->u.rects[i].bptr);
1072       node->u.rects[i].bptr = NULL;
1073     }
1074   /* squeeze the manage flags together */
1075   flag = node->flags.manage & mask;
1076   mask = (~mask) << 1;
1077   node->flags.manage = flag | ((node->flags.manage & mask) >> 1);
1078   /* remove the entry */
1079   for (; i < M_SIZE; i++)
1080     {
1081       node->u.rects[i] = node->u.rects[i + 1];
1082       if (!node->u.rects[i].bptr)
1083         break;
1084     }
1085   if (!node->u.rects[0].bptr)
1086     {
1087       if (node->parent)
1088         __r_delete (node->parent, &node->box);
1089       return true;
1090     }
1091   else
1092     /* propagate boundary adjustment upward */
1093     while (node)
1094       {
1095         adjust_bounds (node);
1096         node = node->parent;
1097       }
1098   return true;
1099 }
1100
1101 bool
1102 r_delete_entry (rtree_t * rtree, const BoxType * box)
1103 {
1104   bool r;
1105
1106   assert (box);
1107   assert (rtree);
1108   r = __r_delete (rtree->root, box);
1109   if (r)
1110     rtree->size--;
1111 #ifdef SLOW_ASSERTS
1112   assert (__r_tree_is_good (rtree->root));
1113 #endif
1114   return r;
1115 }