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