src/rtree.c

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