src/rtree.c

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