Change rtree penalty function in rtree.c to use doubles, not long long.
[geda-pcb/gde.git] / src / rtree.c
blobaf6c8cce72a0184f9e7e1487cd5088c0d1df6937
1 /* $Id$ */
3 /*
4 * COPYRIGHT
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
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.
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.
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.
24 * Contact addresses for paper mail and Email:
25 * harry eaton, 6697 Buttonhole Ct, Columbia, MD 21044 USA
26 * haceaton@aplcomm.jhuapl.edu
30 /* this file, rtree.c, was written and is
31 * Copyright (c) 2004, harry eaton
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.
39 #ifdef HAVE_CONFIG_H
40 #include "config.h"
41 #endif
43 #include "global.h"
45 #include <assert.h>
46 #include <setjmp.h>
48 /* #define DRAWBOX */
50 #ifdef DRAWBOX
51 #include "clip.h"
52 #include "data.h"
53 #endif
54 #include "mymem.h"
56 #include "rtree.h"
58 #ifdef HAVE_LIBDMALLOC
59 #include <dmalloc.h>
60 #endif
62 RCSID ("$Id$");
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.
70 /* the number of entries in each rtree node
71 * 4 - 7 seem to be pretty good settings
73 #define M_SIZE 6
75 /* it seems that sorting the leaf order slows us down
76 * but sometimes gives better routes
78 #undef SORT
79 #define SORT_NONLEAF
81 #define DELETE_BY_POINTER
83 typedef struct
85 const BoxType *bptr; /* pointer to the box */
86 BoxType bounds; /* copy of the box for locality of reference */
87 } Rentry;
89 struct rtree_node
91 BoxType box; /* bounds rectangle of this node */
92 struct rtree_node *parent; /* parent of this node, NULL = root */
93 struct
95 unsigned is_leaf:1; /* this is a leaf node */
96 unsigned manage:31; /* true==should free 'rect.bptr' if node is destroyed */
98 flags;
99 union
101 struct rtree_node *kids[M_SIZE + 1]; /* when not leaf */
102 Rentry rects[M_SIZE + 1]; /* when leaf */
103 } u;
106 #ifdef SLOW_ASSERTS
107 static int
108 __r_node_is_good (struct rtree_node *node)
110 int i, flag = 1;
111 int kind = -1;
112 Boolean last = False;
114 if (node == NULL)
115 return 1;
116 for (i = 0; i < M_SIZE; i++)
118 if (node->flags.is_leaf)
120 if (!node->u.rects[i].bptr)
122 last = True;
123 continue;
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);
152 else
154 if (!node->u.kids[i])
156 last = True;
157 continue;
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);
177 flag <<= 1;
179 /* check that we're completely in the parent's bounds */
180 if (node->parent)
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);
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;
200 /* check the whole tree from this node down for consistency */
201 static Boolean
202 __r_tree_is_good (struct rtree_node *node)
204 int i;
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++)
214 if (!__r_tree_is_good (node->u.kids[i]))
215 return 0;
217 return 1;
219 #endif
220 #ifndef NDEBUG
221 /* print out the tree */
222 void
223 __r_dump_tree (struct rtree_node *node, int depth)
225 int i, j;
226 static int count;
227 static double area;
229 if (depth == 0)
231 area = 0;
232 count = 0;
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)
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);
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
262 else
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++)
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);
316 return;
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);
324 #endif
326 /* Sort the children or entries of a node
327 * according to the largest side.
329 #ifdef SORT
330 static int
331 cmp_box (const BoxType * a, const BoxType * b)
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
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;
355 static void
356 sort_node (struct rtree_node *node)
358 if (node->flags.is_leaf)
360 register Rentry *r, *i, temp;
362 for (r = &node->u.rects[1]; r->bptr; r++)
364 temp = *r;
365 i = r - 1;
366 while (i >= &node->u.rects[0])
368 if (cmp_box (&i->bounds, &r->bounds))
369 break;
370 *(i + 1) = *i;
371 i--;
373 *(i + 1) = temp;
376 #ifdef SORT_NONLEAF
377 else
379 register struct rtree_node **r, **i, *temp;
381 for (r = &node->u.kids[1]; *r; r++)
383 temp = *r;
384 i = r - 1;
385 while (i >= &node->u.kids[0])
387 if (cmp_box (&(*i)->box, &(*r)->box))
388 break;
389 *(i + 1) = *i;
390 i--;
392 *(i + 1) = temp;
395 #endif
397 #else
398 #define sort_node(x)
399 #endif
401 /* set the node bounds large enough to encompass all
402 * of the children's rectangles
404 static void
405 adjust_bounds (struct rtree_node *node)
407 int i;
409 assert (node);
410 assert (node->u.kids[0]);
411 if (node->flags.is_leaf)
413 node->box = node->u.rects[0].bounds;
414 for (i = 1; i < M_SIZE + 1; i++)
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);
424 else
426 node->box = node->u.kids[0]->box;
427 for (i = 1; i < M_SIZE + 1; i++)
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);
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.
444 rtree_t *
445 r_create_tree (const BoxType * boxlist[], int N, int manage)
447 rtree_t *rtree;
448 struct rtree_node *node;
449 int i;
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++)
461 assert (boxlist[i]);
462 r_insert_entry (rtree, boxlist[i], manage);
464 #ifdef SLOW_ASSERTS
465 assert (__r_tree_is_good (rtree->root));
466 #endif
467 return rtree;
470 static void
471 __r_destroy_tree (struct rtree_node *node)
473 int i, flag = 1;
475 if (node->flags.is_leaf)
476 for (i = 0; i < M_SIZE; i++)
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;
484 else
485 for (i = 0; i < M_SIZE; i++)
487 if (!node->u.kids[i])
488 break;
489 __r_destroy_tree (node->u.kids[i]);
491 free (node);
494 /* free the memory associated with an rtree. */
495 void
496 r_destroy_tree (rtree_t ** rtree)
499 __r_destroy_tree ((*rtree)->root);
500 free (*rtree);
501 *rtree = NULL;
504 typedef struct
506 int (*check_it) (const BoxType * region, void *cl);
507 int (*found_it) (const BoxType * box, void *cl);
508 void *closure;
509 } r_arg;
511 /* most of the auto-routing time is spent in this routine
512 * so some careful thought has been given to maximizing the speed
516 __r_search (struct rtree_node *node, const BoxType * query, r_arg * arg)
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->X1);
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.
531 if (node->flags.is_leaf)
533 register int i;
534 if (arg->found_it) /* test this once outside of loop */
536 register int seen = 0;
537 for (i = 0; node->u.rects[i].bptr; i++)
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++;
546 return seen;
548 else
550 register int seen = 0;
551 for (i = 0; node->u.rects[i].bptr; i++)
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++;
559 return seen;
563 /* not a leaf, recurse on lower nodes */
564 if (arg->check_it)
566 int seen = 0;
567 struct rtree_node **n;
568 for (n = &node->u.kids[0]; *n; n++)
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);
577 return seen;
579 else
581 int seen = 0;
582 struct rtree_node **n;
583 for (n = &node->u.kids[0]; *n; n++)
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);
591 return seen;
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
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)
606 r_arg arg;
608 if (!rtree || rtree->size < 1)
609 return 0;
610 if (query)
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);
627 else
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);
636 /*------ r_region_is_empty ------*/
637 static int
638 __r_region_is_empty_rect_in_reg (const BoxType * box, void *cl)
640 jmp_buf *envp = (jmp_buf *) cl;
641 longjmp (*envp, 1); /* found one! */
644 /* return 0 if there are any rectangles in the given region. */
646 r_region_is_empty (rtree_t * rtree, const BoxType * region)
648 jmp_buf env;
649 int r;
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 */
658 struct centroid
660 float x, y, area;
663 /* split the node into two nodes putting clusters in each
664 * use the k-means clustering algorithm
666 struct rtree_node *
667 find_clusters (struct rtree_node *node)
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;
679 for (i = 0; i < M_SIZE + 1; i++)
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);
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)
699 b_X += 10000;
700 a_Y -= 10000;
702 for (tries = 0; tries < M_SIZE; tries++)
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++)
711 float dist1, dist2;
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))
717 belong[i] = True;
718 clust_a++;
720 else
722 belong[i] = False;
723 clust_b++;
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++)
736 if (belong[i])
738 a_X += center[i].x * center[i].area;
739 a_Y += center[i].y * center[i].area;
740 total_a += center[i].area;
742 else
744 b_X += center[i].x * center[i].area;
745 b_Y += center[i].y * center[i].area;
746 total_b += center[i].area;
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;
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)
764 int flag, a_flag, b_flag;
765 flag = a_flag = b_flag = 1;
766 for (i = 0; i < M_SIZE + 1; i++)
768 if (belong[i])
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;
775 else
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;
782 flag <<= 1;
785 else
787 for (i = 0; i < M_SIZE + 1; i++)
789 if (belong[i])
790 node->u.kids[clust_a++] = node->u.kids[i];
791 else
793 node->u.kids[i]->parent = new_node;
794 new_node->u.kids[clust_b++] = node->u.kids[i];
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);
815 /* split a node according to clusters
817 static void
818 split_node (struct rtree_node *node)
820 int i;
821 struct rtree_node *new_node;
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 */
829 struct rtree_node *second;
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;
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)
861 #ifdef SLOW_ASSERTS
862 assert (__r_node_is_good (node->parent));
863 #endif
864 sort_node (node->parent);
865 return;
867 split_node (node->parent);
870 static inline int
871 contained (struct rtree_node *node, const BoxType * query)
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;
880 static inline double
881 penalty (struct rtree_node *node, const BoxType * query)
883 double score;
885 /* Compute the area penalty for inserting here and return.
886 * The penalty is the increase in area necessary
887 * to include the query->
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;
897 static void
898 __r_insert_node (struct rtree_node *node, const BoxType * query,
899 int manage, Boolean force)
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
910 if (node->flags.is_leaf)
912 register int i;
914 if (manage)
916 register int flag = 1;
918 for (i = 0; i < M_SIZE; i++)
920 if (!node->u.rects[i].bptr)
921 break;
922 flag <<= 1;
924 node->flags.manage |= flag;
926 else
928 for (i = 0; i < M_SIZE; i++)
929 if (!node->u.rects[i].bptr)
930 break;
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)
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);
945 if (i < M_SIZE)
947 sort_node (node);
948 return;
950 /* we must split the node */
951 split_node (node);
952 return;
954 else
956 int i;
957 struct rtree_node *best_node;
958 double score, best_score;
960 if (force)
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);
968 /* this node encloses it, but it's not a leaf, so descend the tree */
970 /* First check if any children actually encloses it */
971 assert (node->u.kids[0]);
972 for (i = 0; i < M_SIZE; i++)
974 if (!node->u.kids[i])
975 break;
976 if (contained (node->u.kids[i], query))
978 __r_insert_node (node->u.kids[i], query, manage, False);
979 sort_node (node);
980 return;
984 /* see if there is room for a new leaf node */
985 if (node->u.kids[0]->flags.is_leaf && i < M_SIZE)
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;
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++)
1006 if (!node->u.kids[i])
1007 break;
1008 score = penalty (node->u.kids[i], query);
1009 if (score < best_score)
1011 best_score = score;
1012 best_node = node->u.kids[i];
1015 __r_insert_node (best_node, query, manage, True);
1016 sort_node (node);
1017 return;
1021 void
1022 r_insert_entry (rtree_t * rtree, const BoxType * which, int man)
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++;
1037 Boolean
1038 __r_delete (rtree_t * seed, struct rtree_node *node, const BoxType * query)
1040 int i, flag, mask, a;
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)
1048 for (i = 0; i < M_SIZE; i++)
1050 /* if this is us being removed, free and copy over */
1051 if (node->u.kids[i] == (struct rtree_node *) query)
1053 free ((void *) query);
1054 for (; i < M_SIZE; i++)
1056 node->u.kids[i] = node->u.kids[i + 1];
1057 if (!node->u.kids[i])
1058 break;
1060 /* nobody home here now ? */
1061 if (!node->u.kids[0])
1063 if (!node->parent)
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;
1072 return (__r_delete (seed, node->parent, &node->box));
1074 else
1075 /* propegate boundary adjust upward */
1076 while (node)
1078 adjust_bounds (node);
1079 node = node->parent;
1081 return True;
1083 if (node->u.kids[i])
1085 if (__r_delete (seed, node->u.kids[i], query))
1086 return True;
1088 else
1089 break;
1091 return False;
1093 /* leaf node here */
1094 mask = 0;
1095 a = 1;
1096 for (i = 0; i < M_SIZE; i++)
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;
1110 if (!node->u.rects[i].bptr)
1111 return False; /* not at this leaf */
1112 if (node->flags.manage & a)
1114 free ((void *) node->u.rects[i].bptr);
1115 node->u.rects[i].bptr = NULL;
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++)
1124 node->u.rects[i] = node->u.rects[i + 1];
1125 if (!node->u.rects[i].bptr)
1126 break;
1128 if (!node->u.rects[0].bptr)
1130 if (node->parent)
1131 __r_delete (seed, node->parent, &node->box);
1132 return True;
1134 else
1135 /* propagate boundary adjustment upward */
1136 while (node)
1138 adjust_bounds (node);
1139 node = node->parent;
1141 return True;
1144 Boolean
1145 r_delete_entry (rtree_t * rtree, const BoxType * box)
1147 Boolean r;
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;
1161 __r_sub (struct rtree_node *node, const BoxType * b, const BoxType * a,
1162 jmp_buf * e)
1164 int i;
1166 if (node->flags.is_leaf)
1168 for (i = 0; i < M_SIZE; i++)
1169 if (node->u.rects[i].bptr == b)
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);
1178 return 0;
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;
1187 void
1188 r_substitute (rtree_t * rtree, const BoxType * before, const BoxType * after)
1190 jmp_buf env;
1192 if (before == after)
1193 return;
1194 if (setjmp (env) == 0)
1195 __r_sub (rtree->root, before, after, &env);