n-dimensional rtree

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@51799 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/vector/rtree/index.c

@@ -22,7 +22,7 @@
 #include <assert.h>
 #include <grass/gis.h>
 #include "index.h"
-#include "card.h"
+//#include "card.h"
 
 /* 
  * Make a new index, empty.
@@ -35,9 +35,8 @@ struct RTree *RTreeNewIndex(int fd, off_t rootpos, int ndims)
 {
     struct RTree *new_rtree;
     struct RTree_Node *n;
-    int i;
+    int i, j;
     
-    assert(ndims > 0);
     new_rtree = (struct RTree *)malloc(sizeof(struct RTree));
 
     new_rtree->fd = fd;
@@ -46,19 +45,24 @@ struct RTree *RTreeNewIndex(int fd, off_t rootpos, int ndims)
     new_rtree->nsides = 2 * ndims;
     /* hack to keep compatibility */
     if (ndims < 3)
-	new_rtree->nsides_alloc = 6;
+	new_rtree->ndims_alloc = 3;
     else
-	new_rtree->nsides_alloc = 2 * ndims;
-    
-    new_rtree->rectsize = sizeof(struct RTree_Rect);
+	new_rtree->ndims_alloc = ndims;
+
+    new_rtree->nsides_alloc = 2 * new_rtree->ndims_alloc;
 
     /* init free nodes */
     new_rtree->free_nodes.avail = 0;
     new_rtree->free_nodes.alloc = 0;
     new_rtree->free_nodes.pos = NULL;
 
-    new_rtree->nodesize = sizeof(struct RTree_Node);
-    new_rtree->branchsize = sizeof(struct RTree_Branch);
+    new_rtree->nodesize = sizeof(struct RTree_Node) -
+                          MAXCARD * sizeof(RectReal *) +
+			  MAXCARD * new_rtree->nsides_alloc * sizeof(RectReal);
+
+    new_rtree->branchsize = sizeof(struct RTree_Branch) -
+                            sizeof(RectReal *) +
+			    new_rtree->nsides_alloc * sizeof(RectReal);
     
     /* create empty root node */
     n = RTreeNewNode(new_rtree, 0);
@@ -82,6 +86,15 @@ struct RTree *RTreeNewIndex(int fd, off_t rootpos, int ndims)
 	    new_rtree->used[i][0] = 2;
 	    new_rtree->used[i][1] = 1;
 	    new_rtree->used[i][2] = 0;
+
+	    /* alloc memory for rectangles */
+	    for (j = 0; j < MAXCARD; j++) {
+		RTreeNewRect(&(new_rtree->nb[i][0].n.branch[j].rect), new_rtree);
+		RTreeNewRect(&(new_rtree->nb[i][1].n.branch[j].rect), new_rtree);
+		RTreeNewRect(&(new_rtree->nb[i][2].n.branch[j].rect), new_rtree);
+
+		RTreeNewRect(&(new_rtree->fs[i].sn.branch[j].rect), new_rtree);
+	    }
 	}
 
 	/* write empty root node */
@@ -97,7 +110,6 @@ struct RTree *RTreeNewIndex(int fd, off_t rootpos, int ndims)
 	new_rtree->delete_rect = RTreeDeleteRectF;
 	new_rtree->search_rect = RTreeSearchF;
 	new_rtree->valid_child = RTreeValidChildF;
-	
     }
     else {    /* memory based */
 	new_rtree->nodecard = MAXCARD;
@@ -123,11 +135,29 @@ struct RTree *RTreeNewIndex(int fd, off_t rootpos, int ndims)
     new_rtree->n_nodes = 1;
     new_rtree->n_leafs = 0;
 
+    /* initialize temp variables */
+    RTreeNewRect(&(new_rtree->p.cover[0]), new_rtree);
+    RTreeNewRect(&(new_rtree->p.cover[1]), new_rtree);
+    
+    RTreeNewRect(&(new_rtree->tmpb1.rect), new_rtree);
+    RTreeNewRect(&(new_rtree->tmpb2.rect), new_rtree);
+    RTreeNewRect(&(new_rtree->c.rect), new_rtree);
+    for (i = 0; i <= MAXCARD; i++) {
+	RTreeNewRect(&(new_rtree->BranchBuf[i].rect), new_rtree);
+    }
+    RTreeNewRect(&(new_rtree->rect_0), new_rtree);
+    RTreeNewRect(&(new_rtree->rect_1), new_rtree);
+    RTreeNewRect(&(new_rtree->upperrect), new_rtree);
+    RTreeNewRect(&(new_rtree->orect), new_rtree);
+    new_rtree->center_n = (RectReal *)malloc(new_rtree->ndims_alloc * sizeof(RectReal));
+
     return new_rtree;
 }
 
 void RTreeFreeIndex(struct RTree *t)
 {
+    int i, j;
+
     assert(t);
 
     if (t->fd > -1) {
@@ -137,6 +167,37 @@ void RTreeFreeIndex(struct RTree *t)
     else if (t->root)
 	RTreeDestroyNode(t->root, t->root->level ? t->nodecard : t->leafcard);
 
+    if (t->fd > -1) {  /* file based */
+	/* free node buffer */
+	for (i = 0; i < MAXLEVEL; i++) {
+
+	    /* free memory for rectangles */
+	    for (j = 0; j < MAXCARD; j++) {
+		free(t->nb[i][0].n.branch[j].rect.boundary);
+		free(t->nb[i][1].n.branch[j].rect.boundary);
+		free(t->nb[i][2].n.branch[j].rect.boundary);
+
+		free(t->fs[i].sn.branch[j].rect.boundary);
+	    }
+	}
+    }
+
+    /* free temp variables */
+    free(t->p.cover[0].boundary);
+    free(t->p.cover[1].boundary);
+    
+    free(t->tmpb1.rect.boundary);
+    free(t->tmpb2.rect.boundary);
+    free(t->c.rect.boundary);
+    for (i = 0; i <= MAXCARD; i++) {
+	free(t->BranchBuf[i].rect.boundary);
+    }
+    free(t->rect_0.boundary);
+    free(t->rect_1.boundary);
+    free(t->upperrect.boundary);
+    free(t->orect.boundary);
+    free(t->center_n);
+
     free(t);
     
     return;
@@ -144,7 +205,7 @@ void RTreeFreeIndex(struct RTree *t)
 
 /*
  * Search in an index tree for all data retangles that
- * overlap the argument rectangle.
+ * overlap or touch the argument rectangle.
  * Return the number of qualifying data rects.
  */
 int RTreeSearch(struct RTree *t, struct RTree_Rect *r, SearchHitCallback *shcb,
@@ -228,6 +289,7 @@ void RTreeReInsertNode(struct RTree_Node *n, struct RTree_ListNode **ee)
  */
 void RTreeFreeListBranch(struct RTree_ListBranch *p)
 {
+    free(p->b.rect.boundary);
     free(p);
 }
 

lib/vector/rtree/index.h

@@ -50,6 +50,8 @@ struct RTree_ListBranch
     int level;
 };
 
+/* functions */
+
 /* index.c */
 struct RTree_ListNode *RTreeNewListNode(void);
 void RTreeFreeListNode(struct RTree_ListNode *);
@@ -68,28 +70,36 @@ int RTreeSearchF(struct RTree *, struct RTree_Rect *, SearchHitCallback *,
 		       void *);
 int RTreeInsertRectF(struct RTree_Rect *, union RTree_Child, int, struct RTree *);
 int RTreeDeleteRectF(struct RTree_Rect *, union RTree_Child, struct RTree *);
-int RTreeValidChildF(union RTree_Child *child);
+int RTreeValidChildF(union RTree_Child *);
 
 /* node.c */
-struct RTree_Rect RTreeNodeCover(struct RTree_Node *, struct RTree *);
+void RTreeNodeCover(struct RTree_Node *, struct RTree_Rect *, struct RTree *);
 int RTreeAddBranch(struct RTree_Branch *, struct RTree_Node *, struct RTree_Node **, 
             struct RTree_ListBranch **, struct RTree_Rect *, int *, struct RTree *);
 int RTreePickBranch(struct RTree_Rect *, struct RTree_Node *, struct RTree *);
 void RTreeDisconnectBranch(struct RTree_Node *, int, struct RTree *);
 void RTreePrintNode(struct RTree_Node *, int, struct RTree *);
 void RTreeTabIn(int);
+void RTreeCopyNode(struct RTree_Node *, struct RTree_Node *, struct RTree *);
+void RTreeCopyBranch(struct RTree_Branch *, struct RTree_Branch *, struct RTree *);
 
 /* rect.c */
-void RTreeNewRect(struct RTree_Rect *, struct RTree *);
 void RTreeInitRect(struct RTree_Rect *, struct RTree *);
-struct RTree_Rect RTreeNullRect(void);
+void RTreeNullRect(struct RTree_Rect *, struct RTree *);
 RectReal RTreeRectArea(struct RTree_Rect *, struct RTree *);
 RectReal RTreeRectSphericalVolume(struct RTree_Rect *, struct RTree *);
 RectReal RTreeRectVolume(struct RTree_Rect *, struct RTree *);
 RectReal RTreeRectMargin(struct RTree_Rect *, struct RTree *);
-struct RTree_Rect RTreeCombineRect(struct RTree_Rect *, struct RTree_Rect *, struct RTree *);
+void RTreeCombineRect(struct RTree_Rect *, struct RTree_Rect *, struct RTree_Rect *, struct RTree *);
+void RTreeExpandRect(struct RTree_Rect *, struct RTree_Rect *, struct RTree *);
 int RTreeCompareRect(struct RTree_Rect *, struct RTree_Rect *, struct RTree *);
-void RTreePrintRect(struct RTree_Rect *, int);
+void RTreePrintRect(struct RTree_Rect *, int, struct RTree *);
+
+/*-----------------------------------------------------------------------------
+| Copy second rectangle to first rectangle.
+-----------------------------------------------------------------------------*/
+#define RTreeCopyRect(r1, r2, t) memcpy((r1)->boundary, (r2)->boundary, (t)->nsides_alloc * sizeof(RectReal))
+
 
 /* split.c */
 void RTreeSplitNode(struct RTree_Node *, struct RTree_Branch *, struct RTree_Node *, struct RTree *);

lib/vector/rtree/indexf.c

@@ -24,15 +24,7 @@
 #include <grass/config.h>
 #include <grass/gis.h>
 #include "index.h"
-#include "card.h"
-
-/* stack used for non-recursive insertion/deletion */
-struct fstack
-{
-    struct RTree_Node sn;	/* stack node */
-    int branch_id;	/* branch no to follow down */
-    off_t pos;		/* file position of stack node */
-};
+//#include "card.h"
 
 int RTreeValidChildF(union RTree_Child *child)
 {
@@ -50,8 +42,8 @@ int RTreeSearchF(struct RTree *t, struct RTree_Rect *r,
     struct RTree_Node *n;
     int hitCount = 0, found, currlevel;
     int i;
-    struct fstack s[MAXLEVEL];
     int top = 0;
+    struct fstack *s = t->fs;
 
     /* stack size of t->rootlevel + 1 is enough because of depth first search */
     /* only one node per level on stack at any given time */
@@ -121,12 +113,24 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
 			     struct RTree_ListBranch **ee, int *overflow)
 {
     int i, currlevel;
-    struct RTree_Branch b;
-    struct RTree_Rect nr, *cover;
-    struct RTree_Node *n, *n2, nn;
-    struct fstack s[MAXLEVEL];
+    struct RTree_Rect *cover;
+    struct RTree_Node *n, *n2;
     int top = 0, down = 0;
+    struct RTree_Branch *b = &(t->tmpb2);
     int result;
+    struct fstack *s = t->fs;
+
+    static struct RTree_Node nn;
+    static struct RTree_Rect nr;
+    static int rect_init = 0;
+    
+    if (!rect_init) {
+	for (i = 0; i < MAXCARD; i++) {
+	    RTreeNewRect(&(nn.branch[i].rect), t);
+	}
+	RTreeNewRect(&nr, t);
+	rect_init = t->ndims_alloc;
+    }
 
     n2 = &nn;
 
@@ -149,15 +153,15 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
     
     /* Have reached level for insertion. Add rect, split if necessary */
     if (s[top].sn.level == level) {
-	b.rect = *r;
+	RTreeCopyRect(&(b->rect), r, t);
 	/* child field of leaves contains tid of data record */
-	b.child = child;
+	b->child = child;
 	/* add branch, may split node or remove branches */
 	if (top)
 	    cover = &(s[top - 1].sn.branch[s[top - 1].branch_id].rect);
 	else 
 	    cover = NULL;
-	result = RTreeAddBranch(&b, &(s[top].sn), &n2, ee, cover, overflow, t);
+	result = RTreeAddBranch(b, &(s[top].sn), &n2, ee, cover, overflow, t);
 	/* write out new node if node was split */
 	if (result == 1) {
 	    *newnode_pos = RTreeGetNodePos(t);
@@ -178,7 +182,7 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
 	down = top--;
 	i = s[top].branch_id;
 	if (result == 0) {        /* branch was added */
-	    nr = RTreeCombineRect(r, &(s[top].sn.branch[i].rect), t);
+	    RTreeCombineRect(&(s[top].sn.branch[i].rect), &nr, r, t);
 	    /* rewrite rect */
 	    if (!RTreeCompareRect(&nr, &(s[top].sn.branch[i].rect), t)) {
 		RTreeUpdateRect(&nr, &(s[top].sn), s[top].pos, i, t);
@@ -186,7 +190,7 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
 	}
 	else if (result == 2) {	/* branches were removed */
 	    /* get node cover of previous node */
-	    nr = RTreeNodeCover(&(s[down].sn), t);
+	    RTreeNodeCover(&(s[down].sn), &nr, t);
 	    /* rewrite rect */
 	    if (!RTreeCompareRect(&nr, &(s[top].sn.branch[i].rect), t)) {
 		RTreeUpdateRect(&nr, &(s[top].sn), s[top].pos, i, t);
@@ -194,17 +198,17 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
 	}
 	else if (result == 1) {                /* node was split */
 	    /* get node cover of previous node */
-	    s[top].sn.branch[i].rect = RTreeNodeCover(&(s[down].sn), t);
+	    RTreeNodeCover(&(s[down].sn), &(s[top].sn.branch[i].rect), t);
 	    /* add new branch for new node previously added by RTreeAddBranch() */
-	    b.child.pos = *newnode_pos;
-	    b.rect = RTreeNodeCover(n2, t);
+	    b->child.pos = *newnode_pos;
+	    RTreeNodeCover(n2, &(b->rect), t);
 	    
 	    /* add branch, may split node or remove branches */
 	    if (top)
 		cover = &(s[top - 1].sn.branch[s[top - 1].branch_id].rect);
 	    else
 		cover = NULL;
-	    result = RTreeAddBranch(&b, &(s[top].sn), &n2, ee, cover, overflow, t);
+	    result = RTreeAddBranch(b, &(s[top].sn), &n2, ee, cover, overflow, t);
 	    
 	    /* write out new node if node was split */
 	    if (result == 1) {
@@ -217,7 +221,8 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
 	}
     }
 
-    *newnode = *n2;
+    /* copy node */
+    RTreeCopyNode(newnode, n2, t);
 
     return result;
 }
@@ -233,13 +238,24 @@ static int RTreeInsertRect2F(struct RTree_Rect *r, union RTree_Child child, int
 int RTreeInsertRectF(struct RTree_Rect *r, union RTree_Child child, int level,
                      struct RTree *t)
 {
-    struct RTree_Node oldroot, newroot, newnode;
-    struct RTree_Branch b;
     struct RTree_ListBranch *e, *reInsertList = NULL;
     int result;
     int i, overflow[MAXLEVEL];
     off_t newnode_pos = -1;
 
+    static struct RTree_Node oldroot, newroot, newnode;
+    struct RTree_Branch *b = &(t->tmpb1);
+    static int rect_init = 0;
+    
+    if (!rect_init) {
+	for (i = 0; i < MAXCARD; i++) {
+	    RTreeNewRect(&(oldroot.branch[i].rect), t);
+	    RTreeNewRect(&(newroot.branch[i].rect), t);
+	    RTreeNewRect(&(newnode.branch[i].rect), t);
+	}
+	rect_init = t->ndims_alloc;
+    }
+
     /* R*-tree forced reinsertion: for each level only once */
     for (i = 0; i < MAXLEVEL; i++)
 	overflow[i] = 1;
@@ -255,13 +271,13 @@ int RTreeInsertRectF(struct RTree_Rect *r, union RTree_Child child, int level,
 	RTreeInitNode(t, &newroot, NODETYPE(t->rootlevel, t->fd));
 	newroot.level = t->rootlevel;
 	/* branch for old root */
-	b.rect = RTreeNodeCover(&oldroot, t);
-	b.child.pos = t->rootpos;
-	RTreeAddBranch(&b, &newroot, NULL, NULL, NULL, NULL, t);
+	RTreeNodeCover(&oldroot, &(b->rect), t);
+	b->child.pos = t->rootpos;
+	RTreeAddBranch(b, &newroot, NULL, NULL, NULL, NULL, t);
 	/* branch for new node created by RTreeInsertRect2F() */
-	b.rect = RTreeNodeCover(&newnode, t);
-	b.child.pos = newnode_pos;  /* offset to new node as returned by RTreeInsertRect2F() */
-	RTreeAddBranch(&b, &newroot, NULL, NULL, NULL, NULL, t);
+	RTreeNodeCover(&newnode, &(b->rect), t);
+	b->child.pos = newnode_pos;  /* offset to new node as returned by RTreeInsertRect2F() */
+	RTreeAddBranch(b, &newroot, NULL, NULL, NULL, NULL, t);
 	/* write new root node */
 	t->rootpos = RTreeGetNodePos(t);
 	RTreeWriteNode(&newroot, t);
@@ -270,14 +286,14 @@ int RTreeInsertRectF(struct RTree_Rect *r, union RTree_Child child, int level,
     else if (result == 2) {	/* branches were removed */
 	while (reInsertList) {
 	    /* get next branch in list */
-	    b = reInsertList->b;
+	    RTreeCopyBranch(b, &(reInsertList->b), t);
 	    level = reInsertList->level;
 	    e = reInsertList;
 	    reInsertList = reInsertList->next;
 	    RTreeFreeListBranch(e);
 	    /* reinsert branches */
 	    result =
-		RTreeInsertRect2F(&(b.rect), b.child, level, &newnode, &newnode_pos, t,
+		RTreeInsertRect2F(&(b->rect), b->child, level, &newnode, &newnode_pos, t,
 				 &reInsertList, overflow);
 
 	    if (result == 1) {	/* root split */
@@ -287,13 +303,13 @@ int RTreeInsertRectF(struct RTree_Rect *r, union RTree_Child child, int level,
 		RTreeInitNode(t, &newroot, NODETYPE(t->rootlevel, t->fd));
 		newroot.level = t->rootlevel;
 		/* branch for old root */
-		b.rect = RTreeNodeCover(&oldroot, t);
-		b.child.pos = t->rootpos;
-		RTreeAddBranch(&b, &newroot, NULL, NULL, NULL, NULL, t);
+		RTreeNodeCover(&oldroot, &(b->rect), t);
+		b->child.pos = t->rootpos;
+		RTreeAddBranch(b, &newroot, NULL, NULL, NULL, NULL, t);
 		/* branch for new node created by RTreeInsertRect2F() */
-		b.rect = RTreeNodeCover(&newnode, t);
-		b.child.pos = newnode_pos; 
-		RTreeAddBranch(&b, &newroot, NULL, NULL, NULL, NULL, t);
+		RTreeNodeCover(&newnode, &(b->rect), t);
+		b->child.pos = newnode_pos; 
+		RTreeAddBranch(b, &newroot, NULL, NULL, NULL, NULL, t);
 		/* write new root node */
 		t->rootpos = RTreeGetNodePos(t);
 		RTreeWriteNode(&newroot, t);
@@ -317,13 +333,20 @@ RTreeDeleteRect2F(struct RTree_Rect *r, union RTree_Child child, struct RTree *t
 {
     int i, notfound = 1, currlevel;
     struct RTree_Node *n;
-    struct RTree_Rect nr;
-    struct fstack s[MAXLEVEL];
     int top = 0, down = 0;
     int minfill;
+    struct fstack *s = t->fs;
 
+    static struct RTree_Rect nr;
+    static int rect_init = 0;
+    
     assert(ee);
     
+    if (!rect_init) {
+	RTreeNewRect(&(nr), t);
+	rect_init = 1;
+    }
+
     /* add root node position to stack */
     currlevel = t->rootlevel;
     s[top].pos = t->rootpos;
@@ -384,7 +407,7 @@ RTreeDeleteRect2F(struct RTree_Rect *r, union RTree_Child child, struct RTree *t
 	minfill = (s[down].sn.level ? t->min_node_fill : t->min_leaf_fill);
 	if (s[down].sn.count >= minfill) {
 	    /* just update node cover */
-	    nr = RTreeNodeCover(&(s[down].sn), t);
+	    RTreeNodeCover(&(s[down].sn), &nr, t);
 	    /* rewrite rect */
 	    if (!RTreeCompareRect(&nr, &(s[top].sn.branch[i].rect), t)) {
 		RTreeUpdateRect(&nr, &(s[top].sn), s[top].pos, i, t);
@@ -394,7 +417,8 @@ RTreeDeleteRect2F(struct RTree_Rect *r, union RTree_Child child, struct RTree *t
 	    /* not enough entries in child, eliminate child node */
 	    assert(s[top].sn.branch[i].child.pos == s[down].pos);
 	    n = RTreeNewNode(t, s[down].sn.level);
-	    memcpy(n, &(s[down].sn), t->nodesize);
+	    /* copy node */
+	    RTreeCopyNode(n, &(s[down].sn), t);
 	    RTreeAddNodePos(s[top].sn.branch[i].child.pos, s[down].sn.level, t);
 	    RTreeReInsertNode(n, ee);
 	    RTreeDisconnectBranch(&(s[top].sn), i, t);
@@ -417,9 +441,19 @@ RTreeDeleteRect2F(struct RTree_Rect *r, union RTree_Child child, struct RTree *t
 int RTreeDeleteRectF(struct RTree_Rect *r, union RTree_Child child, struct RTree *t)
 {
     int i;
-    struct RTree_Node *n, rn;
+    struct RTree_Node *n;
     struct RTree_ListNode *e, *reInsertList = NULL;
 
+    static struct RTree_Node rn;
+    static int rect_init = 0;
+    
+    if (!rect_init) {
+	for (i = 0; i < MAXCARD; i++) {
+	    RTreeNewRect(&(rn.branch[i].rect), t);
+	}
+	rect_init = 1;
+    }
+
     if (!RTreeDeleteRect2F(r, child, t, &reInsertList)) {
 	/* found and deleted a data item */
 

lib/vector/rtree/indexm.c

@@ -20,14 +20,7 @@
 #include <assert.h>
 #include <grass/gis.h>
 #include "index.h"
-#include "card.h"
-
-/* stack used for non-recursive insertion/deletion */
-struct stack
-{
-    struct RTree_Node *sn;	/* stack node */
-    int branch_id;	/* branch no to follow down */
-};
+//#include "card.h"
 
 int RTreeValidChildM(union RTree_Child *child)
 {
@@ -45,8 +38,8 @@ int RTreeSearchM(struct RTree *t, struct RTree_Rect *r,
     struct RTree_Node *n;
     int hitCount = 0, found;
     int i;
-    struct stack s[MAXLEVEL];
     int top = 0;
+    struct mstack *s = t->ms;
 
     /* stack size of t->rootlevel + 1 is enough because of depth first search */
     /* only one node per level on stack at any given time */
@@ -113,12 +106,12 @@ static int RTreeInsertRect2M(struct RTree_Rect *r, union RTree_Child child, int
 			     struct RTree_ListBranch **ee, int *overflow)
 {
     int i;
-    struct RTree_Branch b;
     struct RTree_Node *n, *n2;
     struct RTree_Rect *cover;
-    struct stack s[MAXLEVEL];
     int top = 0, down = 0;
     int result;
+    struct RTree_Branch *b = &(t->tmpb2);
+    struct mstack *s = t->ms;
 
     /* add root node to stack */
     s[top].sn = t->root;
@@ -134,15 +127,15 @@ static int RTreeInsertRect2M(struct RTree_Rect *r, union RTree_Child child, int
 
     /* Have reached level for insertion. Remove p rectangles or split */
     if (s[top].sn->level == level) {
-	b.rect = *r;
+	RTreeCopyRect(&(b->rect), r, t);
 	/* child field of leaves contains tid of data record */
-	b.child = child;
+	b->child = child;
 	/* add branch, may split node or remove branches */
 	if (top)
 	    cover = &(s[top - 1].sn->branch[s[top - 1].branch_id].rect);
 	else 
 	    cover = NULL;
-	result = RTreeAddBranch(&b, s[top].sn, &n2, ee, cover, overflow, t);
+	result = RTreeAddBranch(b, s[top].sn, &n2, ee, cover, overflow, t);
 	/* update node count */
 	if (result == 1) {
 	    t->n_nodes++;
@@ -159,19 +152,18 @@ static int RTreeInsertRect2M(struct RTree_Rect *r, union RTree_Child child, int
 	down = top--;
 	i = s[top].branch_id;
 	if (result == 0) {	/* branch was added */
-	    s[top].sn->branch[i].rect =
-		RTreeCombineRect(r, &(s[top].sn->branch[i].rect), t);
+	    RTreeExpandRect(&(s[top].sn->branch[i].rect), r, t);
 	}
 	else if (result == 2) {	/* branches were removed */
 	    /* get node cover of previous node */
-	    s[top].sn->branch[i].rect = RTreeNodeCover(s[down].sn, t);
+	    RTreeNodeCover(s[down].sn, &(s[top].sn->branch[i].rect), t);
 	}
 	else if (result == 1) {	/* node was split */
 	    /* get node cover of previous node */
-	    s[top].sn->branch[i].rect = RTreeNodeCover(s[down].sn, t);
+	    RTreeNodeCover(s[down].sn, &(s[top].sn->branch[i].rect), t);
 	    /* add new branch for new node previously added by RTreeAddBranch() */
-	    b.child.ptr = n2;
-	    b.rect = RTreeNodeCover(b.child.ptr, t);
+	    b->child.ptr = n2;
+	    RTreeNodeCover(b->child.ptr, &(b->rect), t);
 
 	    /* add branch, may split node or remove branches */
 	    if (top)
@@ -179,7 +171,7 @@ static int RTreeInsertRect2M(struct RTree_Rect *r, union RTree_Child child, int
 	    else
 		cover = NULL;
 	    result =
-		RTreeAddBranch(&b, s[top].sn, &n2, ee, cover, overflow, t);
+		RTreeAddBranch(b, s[top].sn, &n2, ee, cover, overflow, t);
 
 	    /* update node count */
 	    if (result == 1) {
@@ -205,11 +197,11 @@ int RTreeInsertRectM(struct RTree_Rect *r, union RTree_Child child, int level,
 		     struct RTree *t)
 {
     struct RTree_Node *newnode, *newroot;
-    struct RTree_Branch b;
     struct RTree_ListBranch *reInsertList = NULL;
     struct RTree_ListBranch *e;
     int result;
     int i, overflow[MAXLEVEL];
+    struct RTree_Branch *b = &(t->tmpb1);
 
     /* R*-tree forced reinsertion: for each level only once */
     for (i = 0; i < MAXLEVEL; i++)
@@ -225,13 +217,13 @@ int RTreeInsertRectM(struct RTree_Rect *r, union RTree_Child child, int level,
 	newroot = RTreeNewNode(t, t->rootlevel);
 	newroot->level = t->rootlevel;
 	/* branch for old root */
-	b.rect = RTreeNodeCover(t->root, t);
-	b.child.ptr = t->root;
-	RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+	RTreeNodeCover(t->root, &(b->rect), t);
+	b->child.ptr = t->root;
+	RTreeAddBranch(b, newroot, NULL, NULL, NULL, NULL, t);
 	/* branch for new node created by RTreeInsertRect2() */
-	b.rect = RTreeNodeCover(newnode, t);
-	b.child.ptr = newnode;
-	RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+	RTreeNodeCover(newnode, &(b->rect), t);
+	b->child.ptr = newnode;
+	RTreeAddBranch(b, newroot, NULL, NULL, NULL, NULL, t);
 	/* set new root node */
 	t->root = newroot;
 	t->n_nodes++;
@@ -239,14 +231,14 @@ int RTreeInsertRectM(struct RTree_Rect *r, union RTree_Child child, int level,
     else if (result == 2) {	/* branches were removed */
 	while (reInsertList) {
 	    /* get next branch in list */
-	    b = reInsertList->b;
+	    RTreeCopyBranch(b, &(reInsertList->b), t);
 	    level = reInsertList->level;
 	    e = reInsertList;
 	    reInsertList = reInsertList->next;
 	    RTreeFreeListBranch(e);
 	    /* reinsert branches */
 	    result =
-		RTreeInsertRect2M(&(b.rect), b.child, level, &newnode, t,
+		RTreeInsertRect2M(&(b->rect), b->child, level, &newnode, t,
 				 &reInsertList, overflow);
 
 	    if (result == 1) {	/* root split */
@@ -255,13 +247,13 @@ int RTreeInsertRectM(struct RTree_Rect *r, union RTree_Child child, int level,
 		newroot = RTreeNewNode(t, t->rootlevel);
 		newroot->level = t->rootlevel;
 		/* branch for old root */
-		b.rect = RTreeNodeCover(t->root, t);
-		b.child.ptr = t->root;
-		RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+		RTreeNodeCover(t->root, &(b->rect), t);
+		b->child.ptr = t->root;
+		RTreeAddBranch(b, newroot, NULL, NULL, NULL, NULL, t);
 		/* branch for new node created by RTreeInsertRect2() */
-		b.rect = RTreeNodeCover(newnode, t);
-		b.child.ptr = newnode; 
-		RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+		RTreeNodeCover(newnode, &(b->rect), t);
+		b->child.ptr = newnode; 
+		RTreeAddBranch(b, newroot, NULL, NULL, NULL, NULL, t);
 		/* set new root node */
 		t->root = newroot;
 		t->n_nodes++;
@@ -284,9 +276,9 @@ RTreeDeleteRect2M(struct RTree_Rect *r, union RTree_Child child, struct RTree *t
 {
     int i, notfound = 1;
     struct RTree_Node *n;
-    struct stack s[MAXLEVEL];
     int top = 0, down = 0;
     int minfill;
+    struct mstack *s = t->ms;
 
     assert(ee);
 
@@ -348,7 +340,7 @@ RTreeDeleteRect2M(struct RTree_Rect *r, union RTree_Child child, struct RTree *t
 	minfill = (s[down].sn->level ? t->min_node_fill : t->min_leaf_fill);
 	if (s[down].sn->count >= minfill) {
 	    /* just update node cover */
-	    s[top].sn->branch[i].rect = RTreeNodeCover(s[down].sn, t);
+	    RTreeNodeCover(s[down].sn, &(s[top].sn->branch[i].rect), t);
 	}
 	else {
 	    /* not enough entries in child, eliminate child node */

lib/vector/rtree/io.c

@@ -67,11 +67,32 @@ off_t RTreeGetNodePos(struct RTree *t)
     }
 }
 
+/* read branch from file */
+size_t RTreeReadBranch(struct RTree_Branch *b, struct RTree *t)
+{
+    size_t size = 0;
+
+    size += read(t->fd, b->rect.boundary, t->nsides * sizeof(RectReal));
+    size += read(t->fd, &(b->child), sizeof(union RTree_Child));
+
+    return size;
+}
+
 /* read node from file */
 size_t RTreeReadNode(struct RTree_Node *n, off_t nodepos, struct RTree *t)
 {
+    int i;
+    size_t size = 0;
+
     lseek(t->fd, nodepos, SEEK_SET);
-    return read(t->fd, n, t->nodesize);
+    size += read(t->fd, &(n->count), sizeof(int));
+    size += read(t->fd, &(n->level), sizeof(int));
+
+    for (i = 0; i < MAXCARD; i++) {
+	size += RTreeReadBranch(&(n->branch[i]), t);
+    }
+
+    return size;
 }
 
 /* get node from buffer or file */
@@ -100,14 +121,36 @@ void RTreeGetNode(struct RTree_Node *n, off_t nodepos, int level, struct RTree *
 	t->used[level][1] = t->used[level][0];
 	t->used[level][0] = which; 
     }
-    *n = t->nb[level][which].n;
+    /* copy node */
+    RTreeCopyNode(n, &(t->nb[level][which].n), t);
+}
+
+/* write branch to file */
+size_t RTreeWriteBranch(struct RTree_Branch *b, struct RTree *t)
+{
+    size_t size = 0;
+
+    size += write(t->fd, b->rect.boundary, t->nsides * sizeof(RectReal));
+    size += write(t->fd, &(b->child), sizeof(union RTree_Child));
+
+    return size;
 }
 
 /* write new node to file */
 size_t RTreeWriteNode(struct RTree_Node *n, struct RTree *t)
 {
+    int i;
+    size_t size = 0;
+
     /* file position must be set first with RTreeGetFNodePos() */
-    return write(t->fd, n, t->nodesize);
+    size += write(t->fd, &(n->count), sizeof(int));
+    size += write(t->fd, &(n->level), sizeof(int));
+
+    for (i = 0; i < MAXCARD; i++) {
+	size += RTreeWriteBranch(&(n->branch[i]), t);
+    }
+
+    return size;
 }
 
 /* rewrite updated node to file */
@@ -122,7 +165,8 @@ void RTreePutNode(struct RTree_Node *n, off_t nodepos, struct RTree *t)
 {
     int which = (nodepos == t->nb[n->level][2].pos ? 2 : nodepos == t->nb[n->level][1].pos);
     
-    t->nb[n->level][which].n = *n;
+    /* copy node */
+    RTreeCopyNode(&(t->nb[n->level][which].n), n, t);
     t->nb[n->level][which].dirty = 1;
 
     /* make it mru */
@@ -140,9 +184,15 @@ void RTreePutNode(struct RTree_Node *n, off_t nodepos, struct RTree *t)
 /* update rectangle */
 void RTreeUpdateRect(struct RTree_Rect *r, struct RTree_Node *n, off_t nodepos, int b, struct RTree *t)
 {
-    int which = (nodepos == t->nb[n->level][2].pos ? 2 : nodepos == t->nb[n->level][1].pos);
+    int i;
+    int which = (nodepos == t->nb[n->level][2].pos ?
+		 2 : nodepos == t->nb[n->level][1].pos);
     
-    t->nb[n->level][which].n.branch[b].rect = n->branch[b].rect = *r;
+    for (i = 0; i < t->nsides; i++) {
+	t->nb[n->level][which].n.branch[b].rect.boundary[i] =
+	                  n->branch[b].rect.boundary[i] = r->boundary[i];
+    }
+
     t->nb[n->level][which].dirty = 1;
 
     /* make it mru */

lib/vector/rtree/node.c

@@ -18,9 +18,11 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include <assert.h>
 #include <grass/gis.h>
 #include "index.h"
+#include "split.h"
 #include "card.h"
 
 /* rectangle distances for forced removal */
@@ -31,27 +33,27 @@ struct dist
 };
 
 /* Initialize one branch cell in an internal node. */
-static void RTreeInitNodeBranchM(struct RTree *t, struct RTree_Branch *b)
+static void RTreeInitNodeBranchM(struct RTree_Branch *b, struct RTree *t)
 {
     RTreeInitRect(&(b->rect), t);
     b->child.ptr = NULL;
 }
 
 /* Initialize one branch cell in an internal node. */
-static void RTreeInitNodeBranchF(struct RTree *t, struct RTree_Branch *b)
+static void RTreeInitNodeBranchF(struct RTree_Branch *b, struct RTree *t)
 {
     RTreeInitRect(&(b->rect), t);
     b->child.pos = -1;
 }
 
 /* Initialize one branch cell in a leaf node. */
-static void RTreeInitLeafBranch(struct RTree *t, struct RTree_Branch *b)
+static void RTreeInitLeafBranch(struct RTree_Branch *b, struct RTree *t)
 {
     RTreeInitRect(&(b->rect), t);
     b->child.id = 0;
 }
 
-static void (*RTreeInitBranch[3]) () = {
+static void (*RTreeInitBranch[3]) (struct RTree_Branch *b, struct RTree *t) = {
     RTreeInitLeafBranch, RTreeInitNodeBranchM, RTreeInitNodeBranchF
 };
 
@@ -65,45 +67,79 @@ void RTreeInitNode(struct RTree *t, struct RTree_Node *n, int type)
     n->level = -1;
 
     for (i = 0; i < MAXCARD; i++)
-	RTreeInitBranch[type](t, &(n->branch[i]));
+	RTreeInitBranch[type](&(n->branch[i]), t);
 }
 
 /* Make a new node and initialize to have all branch cells empty. */
 struct RTree_Node *RTreeNewNode(struct RTree *t, int level)
 {
+    int i;
     struct RTree_Node *n;
 
-    n = (struct RTree_Node *)malloc((size_t) t->nodesize);
+    n = (struct RTree_Node *)malloc(sizeof(struct RTree_Node));
     assert(n);
-    RTreeInitNode(t, n, NODETYPE(level, t->fd));
+
+    n->count = 0;
+    n->level = -1;
+
+    for (i = 0; i < MAXCARD; i++) {
+	RTreeNewRect(&(n->branch[i].rect), t);
+	RTreeInitBranch[NODETYPE(level, t->fd)](&(n->branch[i]), t);
+    }
+
     return n;
 }
 
 void RTreeFreeNode(struct RTree_Node *n)
 {
+    int i;
+
     assert(n);
+
+    for (i = 0; i < MAXCARD; i++)
+	free(n->branch[i].rect.boundary);
+
     free(n);
 }
 
+/* copy node 2 to node 1 */
+void RTreeCopyNode(struct RTree_Node *n1, struct RTree_Node *n2, struct RTree *t)
+{
+    int i;
+
+    assert(n1 && n2);
+
+    n1->count = n2->count;
+    n1->level = n2->level;
+    for (i = 0; i < MAXCARD; i++) {
+	RTreeCopyBranch(&(n1->branch[i]), &(n2->branch[i]), t);
+    }
+}
+
+/* copy branch 2 to branch 1 */
+void RTreeCopyBranch(struct RTree_Branch *b1, struct RTree_Branch *b2, struct RTree *t)
+{
+    b1->child = b2->child;
+    RTreeCopyRect(&(b1->rect), &(b2->rect), t);
+}
+
 /*
  * Find the smallest rectangle that includes all rectangles in
  * branches of a node.
  */
-struct RTree_Rect RTreeNodeCover(struct RTree_Node *n, struct RTree *t)
+void RTreeNodeCover(struct RTree_Node *n, struct RTree_Rect *r, struct RTree *t)
 {
     int i, first_time = 1;
-    struct RTree_Rect r;
 
-    RTreeInitRect(&r, t);
     if ((n)->level > 0) { /* internal node */
 	for (i = 0; i < t->nodecard; i++) {
 	    if (t->valid_child(&(n->branch[i].child))) {
 		if (first_time) {
-		    r = n->branch[i].rect;
+		    RTreeCopyRect(r, &(n->branch[i].rect), t);
 		    first_time = 0;
 		}
 		else
-		    r = RTreeCombineRect(&r, &(n->branch[i].rect), t);
+		    RTreeExpandRect(r, &(n->branch[i].rect), t);
 	    }
 	}
     }
@@ -111,15 +147,14 @@ struct RTree_Rect RTreeNodeCover(struct RTree_Node *n, struct RTree *t)
 	for (i = 0; i < t->leafcard; i++) {
 	    if (n->branch[i].child.id) {
 		if (first_time) {
-		    r = n->branch[i].rect;
+		    RTreeCopyRect(r, &(n->branch[i].rect), t);
 		    first_time = 0;
 		}
 		else
-		    r = RTreeCombineRect(&r, &(n->branch[i].rect), t);
+		    RTreeExpandRect(r, &(n->branch[i].rect), t);
 	    }
 	}
     }
-    return r;
 }
 
 /*
@@ -140,7 +175,6 @@ static int RTreePickLeafBranch(struct RTree_Rect *r, struct RTree_Node *n, struc
     int i, j;
     RectReal increase, bestIncr = -1, area, bestArea = 0;
     int best = 0, bestoverlap;
-    struct RTree_Rect tmp_rect;
     int overlap;
 
     bestoverlap = t->nodecard + 1;
@@ -150,15 +184,15 @@ static int RTreePickLeafBranch(struct RTree_Rect *r, struct RTree_Node *n, struc
     for (i = 0; i < t->nodecard; i++) {
 	if (t->valid_child(&(n->branch[i].child))) {
 	    rr = &n->branch[i].rect;
-	    tmp_rect = RTreeCombineRect(r, rr, t);
+	    RTreeCombineRect(r, rr, &(t->orect), t);
 	    area = RTreeRectSphericalVolume(rr, t);
-	    increase = RTreeRectSphericalVolume(&tmp_rect, t) - area;
+	    increase = RTreeRectSphericalVolume(&(t->orect), t) - area;
 
 	    overlap = 0;
 	    for (j = 0; j < t->leafcard; j++) {
 		if (j != i) {
 		    rr = &n->branch[j].rect;
-		    overlap += RTreeOverlap(&tmp_rect, rr, t);
+		    overlap += RTreeOverlap(&(t->orect), rr, t);
 		}
 	    }
 
@@ -199,7 +233,6 @@ int RTreePickBranch(struct RTree_Rect *r, struct RTree_Node *n, struct RTree *t)
     int i, first_time = 1;
     RectReal increase, bestIncr = (RectReal) -1, area, bestArea = 0;
     int best = 0;
-    struct RTree_Rect tmp_rect;
 
     assert((n)->level > 0);	/* must not be called on leaf node */
 
@@ -210,8 +243,8 @@ int RTreePickBranch(struct RTree_Rect *r, struct RTree_Node *n, struct RTree *t)
 	if (t->valid_child(&(n->branch[i].child))) {
 	    rr = &n->branch[i].rect;
 	    area = RTreeRectSphericalVolume(rr, t);
-	    tmp_rect = RTreeCombineRect(r, rr, t);
-	    increase = RTreeRectSphericalVolume(&tmp_rect, t) - area;
+	    RTreeCombineRect(r, rr, &(t->orect), t);
+	    increase = RTreeRectSphericalVolume(&(t->orect), t) - area;
 	    if (increase < bestIncr || first_time) {
 		best = i;
 		bestArea = area;
@@ -234,14 +267,14 @@ void RTreeDisconnectBranch(struct RTree_Node *n, int i, struct RTree *t)
 	assert(n && i >= 0 && i < t->nodecard);
 	assert(t->valid_child(&(n->branch[i].child)));
 	if (t->fd < 0)
-	    RTreeInitNodeBranchM(t, &(n->branch[i]));
+	    RTreeInitNodeBranchM(&(n->branch[i]), t);
 	else
-	    RTreeInitNodeBranchF(t, &(n->branch[i]));
+	    RTreeInitNodeBranchF(&(n->branch[i]), t);
     }
     else {
 	assert(n && i >= 0 && i < t->leafcard);
 	assert(n->branch[i].child.id);
-	RTreeInitLeafBranch(t, &(n->branch[i]));
+	RTreeInitLeafBranch(&(n->branch[i]), t);
     }
 
     n->count--;
@@ -398,9 +431,15 @@ static void RTreeQuicksortDist(struct dist *d, int n)
  * Allocate space for a branch in the list used in InsertRect to
  * store branches of nodes that are too full.
  */
-static struct RTree_ListBranch *RTreeNewListBranch(void)
+static struct RTree_ListBranch *RTreeNewListBranch(struct RTree *t)
 {
-    return (struct RTree_ListBranch *)malloc(sizeof(struct RTree_ListBranch));
+    struct RTree_ListBranch *p = 
+       (struct RTree_ListBranch *)malloc(sizeof(struct RTree_ListBranch));
+
+    assert(p);
+    RTreeNewRect(&(p->b.rect), t);
+
+    return p;
 }
 
 /* 
@@ -408,12 +447,12 @@ static struct RTree_ListBranch *RTreeNewListBranch(void)
  * be reinserted into the index structure.
  */
 static void RTreeReInsertBranch(struct RTree_Branch b, int level,
-				struct RTree_ListBranch **ee)
+				struct RTree_ListBranch **ee, struct RTree *t)
 {
     register struct RTree_ListBranch *l;
 
-    l = RTreeNewListBranch();
-    l->b = b;
+    l = RTreeNewListBranch(t);
+    RTreeCopyBranch(&(l->b), &b, t);
     l->level = level;
     l->next = *ee;
     *ee = l;
@@ -429,10 +468,11 @@ static void RTreeRemoveBranches(struct RTree_Node *n, struct RTree_Branch *b,
 				struct RTree *t)
 {
     int i, j, maxkids, type;
-    RectReal center_n[NUMDIMS], center_r, delta;
-    struct RTree_Branch branchbuf[MAXCARD + 1];
+    RectReal center_r, delta;
     struct dist rdist[MAXCARD + 1];
-    struct RTree_Rect new_cover;
+
+    struct RTree_Rect *new_cover = &(t->orect);
+    RectReal *center_n = t->center_n;
 
     assert(cover);
 
@@ -441,35 +481,35 @@ static void RTreeRemoveBranches(struct RTree_Node *n, struct RTree_Branch *b,
     
     assert(n->count == maxkids);	/* must be full */
 
-    new_cover = RTreeCombineRect(cover, &(b->rect), t);
+    RTreeCombineRect(cover, &(b->rect), new_cover, t);
 
     /* center coords of node cover */
     for (j = 0; j < t->ndims; j++) {
-	center_n[j] = (new_cover.boundary[j + NUMDIMS] + new_cover.boundary[j]) / 2;
+	center_n[j] = (new_cover->boundary[j + t->ndims_alloc] + new_cover->boundary[j]) / 2;
     }
 
     /* compute distances of child rectangle centers to node cover center */
     for (i = 0; i < maxkids; i++) {
-	branchbuf[i] = n->branch[i];
+	RTreeCopyBranch(&(t->BranchBuf[i]), &(n->branch[i]), t);
 	rdist[i].distance = 0;
 	rdist[i].id = i;
 	for (j = 0; j < t->ndims; j++) {
 	    center_r =
-		(branchbuf[i].rect.boundary[j + NUMDIMS] +
-		 branchbuf[i].rect.boundary[j]) / 2;
+		(t->BranchBuf[i].rect.boundary[j + t->ndims_alloc] +
+		 t->BranchBuf[i].rect.boundary[j]) / 2;
 	    delta = center_n[j] - center_r;
 	    rdist[i].distance += delta * delta;
 	}
 
-	RTreeInitBranch[type](t, &(n->branch[i]));
+	RTreeInitBranch[type](&(n->branch[i]), t);
     }
 
     /* new branch */
-    branchbuf[maxkids] = *b;
+    RTreeCopyBranch(&(t->BranchBuf[maxkids]), b, t);
     rdist[maxkids].distance = 0;
     for (j = 0; j < t->ndims; j++) {
 	center_r =
-	    (b->rect.boundary[j + NUMDIMS] +
+	    (b->rect.boundary[j + t->ndims_alloc] +
 	     b->rect.boundary[j]) / 2;
 	delta = center_n[j] - center_r;
 	rdist[maxkids].distance += delta * delta;
@@ -481,11 +521,11 @@ static void RTreeRemoveBranches(struct RTree_Node *n, struct RTree_Branch *b,
 
     /* put largest three in branch list, farthest from center first */
     for (i = 0; i < FORCECARD; i++) {
-	RTreeReInsertBranch(branchbuf[rdist[maxkids - i].id], n->level, ee);
+	RTreeReInsertBranch(t->BranchBuf[rdist[maxkids - i].id], n->level, ee, t);
     }
     /* put remaining in node, closest to center first */
     for (i = 0; i < maxkids - FORCECARD + 1; i++) {
-	n->branch[i] = branchbuf[rdist[i].id];
+	RTreeCopyBranch(&(n->branch[i]), &(t->BranchBuf[rdist[i].id]), t);
     }
     n->count = maxkids - FORCECARD + 1;
 }
@@ -509,7 +549,9 @@ int RTreeAddBranch(struct RTree_Branch *b, struct RTree_Node *n,
 	if ((n)->level > 0) {   /* internal node */
 	    for (i = 0; i < maxkids; i++) {	/* find empty branch */
 		if (!t->valid_child(&(n->branch[i].child))) {
-		    n->branch[i] = *b;
+		    /* copy branch */
+		    n->branch[i].child = b->child;
+		    RTreeCopyRect(&(n->branch[i].rect), &(b->rect), t);
 		    n->count++;
 		    break;
 		}
@@ -519,7 +561,9 @@ int RTreeAddBranch(struct RTree_Branch *b, struct RTree_Node *n,
 	else if ((n)->level == 0) {   /* leaf */
 	    for (i = 0; i < maxkids; i++) {	/* find empty branch */
 		if (n->branch[i].child.id == 0) {
-		    n->branch[i] = *b;
+		    /* copy branch */
+		    n->branch[i].child = b->child;
+		    RTreeCopyRect(&(n->branch[i].rect), &(b->rect), t);
 		    n->count++;
 		    break;
 		}
@@ -563,7 +607,7 @@ void RTreeTabIn(int depth)
 
 static void RTreePrintBranch(struct RTree_Branch *b, int depth, struct RTree *t)
 {
-    RTreePrintRect(&(b->rect), depth);
+    RTreePrintRect(&(b->rect), depth, t);
     RTreePrintNode(b->child.ptr, depth, t);
 }
 
@@ -588,7 +632,7 @@ void RTreePrintNode(struct RTree_Node *n, int depth, struct RTree *t)
     for (i = 0; i < maxkids; i++) {
 	if (n->level == 0) {
 	    RTreeTabIn(depth);
-	    RTreePrintRect(&(n->branch[i].rect), depth);
+	    RTreePrintRect(&(n->branch[i].rect), depth, t);
 	    fprintf(stdout, "\t%d: data id = %d\n", i,
 		  n->branch[i].child.id);
 	}

lib/vector/rtree/rect.c

@@ -27,23 +27,19 @@
 #define BIG_NUM (FLT_MAX/4.0)
 
 
-#define Undefined(x) ((x)->boundary[0] > (x)->boundary[NUMDIMS])
+#define Undefined(x, t) ((x)->boundary[0] > (x)->boundary[t->ndims_alloc])
 #define MIN(a, b) ((a) < (b) ? (a) : (b))
 #define MAX(a, b) ((a) > (b) ? (a) : (b))
 
 
-#if 0
 /*-----------------------------------------------------------------------------
 | Create a new rectangle for a given tree
 -----------------------------------------------------------------------------*/
 void RTreeNewRect(struct RTree_Rect *r, struct RTree *t)
 {
-    r->boundary = (struct RTree_Rect *)malloc((size_t) t->nsides_alloc);
+    r->boundary = (RectReal *)malloc((size_t) t->nsides_alloc * sizeof(RectReal));
     assert(r->boundary);
-
-    RTreeInitRect(r, t);
 }
-#endif
 
 /*-----------------------------------------------------------------------------
 | Initialize a rectangle to have all 0 coordinates.
@@ -52,25 +48,26 @@ void RTreeInitRect(struct RTree_Rect *r, struct RTree *t)
 {
     register int i;
 
-    for (i = 0; i < t->nsides; i++)
+    for (i = 0; i < t->nsides_alloc; i++)
 	r->boundary[i] = (RectReal) 0;
 }
 
-
 /*-----------------------------------------------------------------------------
 | Return a rect whose first low side is higher than its opposite side -
 | interpreted as an undefined rect.
 -----------------------------------------------------------------------------*/
-struct RTree_Rect RTreeNullRect(void)
+void RTreeNullRect(struct RTree_Rect *r, struct RTree *t)
 {
-    struct RTree_Rect r;
     register int i;
 
-    r.boundary[0] = (RectReal) 1;
-    r.boundary[NUMDIMS] = (RectReal) - 1;
-    for (i = 1; i < NUMDIMS; i++)
-	r.boundary[i] = r.boundary[i + NUMDIMS] = (RectReal) 0;
-    return r;
+    assert(r);
+
+    r->boundary[0] = (RectReal) 1;
+    r->boundary[t->nsides - 1] = (RectReal) - 1;
+    for (i = 1; i < t->ndims; i++)
+	r->boundary[i] = r->boundary[i + t->ndims_alloc] = (RectReal) 0;
+
+    return;
 }
 
 
@@ -139,7 +136,7 @@ void RTreeSearchRect(struct RTree_Rect *Search, struct RTree_Rect *Data)
 /*-----------------------------------------------------------------------------
 | Print out the data for a rectangle.
 -----------------------------------------------------------------------------*/
-void RTreePrintRect(struct RTree_Rect *R, int depth)
+void RTreePrintRect(struct RTree_Rect *R, int depth, struct RTree *t)
 {
     register struct RTree_Rect *r = R;
     register int i;
@@ -148,9 +145,9 @@ void RTreePrintRect(struct RTree_Rect *R, int depth)
 
     RTreeTabIn(depth);
     fprintf(stdout, "rect:\n");
-    for (i = 0; i < NUMDIMS; i++) {
+    for (i = 0; i < t->ndims_alloc; i++) {
 	RTreeTabIn(depth + 1);
-	fprintf(stdout, "%f\t%f\n", r->boundary[i], r->boundary[i + NUMDIMS]);
+	fprintf(stdout, "%f\t%f\n", r->boundary[i], r->boundary[i + t->ndims_alloc]);
     }
 }
 
@@ -164,12 +161,13 @@ RectReal RTreeRectVolume(struct RTree_Rect *R, struct RTree *t)
     register RectReal volume = (RectReal) 1;
 
     assert(r);
-    if (Undefined(r))
+    if (Undefined(r, t))
 	return (RectReal) 0;
 
     for (i = 0; i < t->ndims; i++)
-	volume *= r->boundary[i + NUMDIMS] - r->boundary[i];
+	volume *= r->boundary[i + t->ndims_alloc] - r->boundary[i];
     assert(volume >= 0.0);
+
     return volume;
 }
 
@@ -195,7 +193,7 @@ static double sphere_volume(double dimension)
     log_volume = dimension / 2.0 * log(M_PI) - log_gamma;
     return exp(log_volume);
 }
-static const double UnitSphereVolume = sphere_volume(NUMDIMS);
+static const double UnitSphereVolume = sphere_volume(20);
 
 #else
 
@@ -248,7 +246,7 @@ RectReal RTreeRectSphericalVolume(struct RTree_Rect *R, struct RTree *t)
     RectReal maxsize = (RectReal) 0, c_size;
 
     assert(r);
-    if (Undefined(r))
+    if (Undefined(r, t))
 	return (RectReal) 0;
     for (i = 0; i < t->ndims; i++) {
 	c_size = r->boundary[i + NUMDIMS] - r->boundary[i];
@@ -269,10 +267,10 @@ RectReal RTreeRectSphericalVolume(struct RTree_Rect *r, struct RTree *t)
     double sum_of_squares = 0, radius, half_extent;
 
     assert(r);
-    if (Undefined(r))
+    if (Undefined(r, t))
 	return (RectReal) 0;
     for (i = 0; i < t->ndims; i++) {
-	half_extent = (r->boundary[i + NUMDIMS] - r->boundary[i]) / 2;
+	half_extent = (r->boundary[i + t->ndims_alloc] - r->boundary[i]) / 2;
 
 	sum_of_squares += half_extent * half_extent;
     }
@@ -290,7 +288,7 @@ RectReal RTreeRectSurfaceArea(struct RTree_Rect *r, struct RTree *t)
     RectReal j_extent, sum = (RectReal) 0;
 
     assert(r);
-    if (Undefined(r))
+    if (Undefined(r, t))
 	return (RectReal) 0;
 
     for (i = 0; i < t->ndims; i++) {
@@ -299,7 +297,7 @@ RectReal RTreeRectSurfaceArea(struct RTree_Rect *r, struct RTree *t)
 	for (j = 0; j < t->ndims; j++)
 	    /* exclude i extent from product in this dimension */
 	    if (i != j) {
-		j_extent = r->boundary[j + NUMDIMS] - r->boundary[j];
+		j_extent = r->boundary[j + t->ndims_alloc] - r->boundary[j];
 
 		face_area *= j_extent;
 	    }
@@ -321,7 +319,7 @@ RectReal RTreeRectMargin(struct RTree_Rect *r, struct RTree *t)
     assert(r);
 
     for (i = 0; i < t->ndims; i++) {
-	margin += r->boundary[i + NUMDIMS] - r->boundary[i];
+	margin += r->boundary[i + t->ndims_alloc] - r->boundary[i];
     }
 
     return margin;
@@ -331,33 +329,64 @@ RectReal RTreeRectMargin(struct RTree_Rect *r, struct RTree *t)
 /*-----------------------------------------------------------------------------
 | Combine two rectangles, make one that includes both.
 -----------------------------------------------------------------------------*/
-struct RTree_Rect RTreeCombineRect(struct RTree_Rect *r, struct RTree_Rect *rr, struct RTree *t)
+void RTreeCombineRect(struct RTree_Rect *r1, struct RTree_Rect *r2,
+		      struct RTree_Rect *r3, struct RTree *t)
 {
     int i, j;
-    struct RTree_Rect new_rect;
 
-    assert(r && rr);
+    assert(r1 && r2 && r3);
 
-    if (Undefined(r))
-	return *rr;
+    if (Undefined(r1, t)) {
+	for (i = 0; i < t->nsides; i++)
+	    r3->boundary[i] = r2->boundary[i];
+	
+	return;
+    }
 
-    if (Undefined(rr))
-	return *r;
+    if (Undefined(r2, t)) {
+	for (i = 0; i < t->nsides; i++)
+	    r3->boundary[i] = r1->boundary[i];
+	
+	return;
+    }
 
     for (i = 0; i < t->ndims; i++) {
-	new_rect.boundary[i] = MIN(r->boundary[i], rr->boundary[i]);
-	j = i + NUMDIMS;
-	new_rect.boundary[j] = MAX(r->boundary[j], rr->boundary[j]);
+	r3->boundary[i] = MIN(r1->boundary[i], r2->boundary[i]);
+	j = i + t->ndims_alloc;
+	r3->boundary[j] = MAX(r1->boundary[j], r2->boundary[j]);
+    }
+    for (i = t->ndims; i < t->ndims_alloc; i++) {
+	r3->boundary[i] = 0;
+	j = i + t->ndims_alloc;
+	r3->boundary[j] = 0;
     }
+}
+
 
-    /* set remaining boundaries to zero */
-    for (; i < NUMDIMS; i++) {
-	new_rect.boundary[i] = MIN(r->boundary[i], rr->boundary[i]);
-	j = i + NUMDIMS;
-	new_rect.boundary[j] = MAX(r->boundary[j], rr->boundary[j]);
+/*-----------------------------------------------------------------------------
+| Expand first rectangle to cover second rectangle.
+-----------------------------------------------------------------------------*/
+void RTreeExpandRect(struct RTree_Rect *r1, struct RTree_Rect *r2,
+		     struct RTree *t)
+{
+    int i, j;
+
+    assert(r1 && r2);
+
+    if (Undefined(r2, t))
+	return;
+
+    for (i = 0; i < t->ndims; i++) {
+	r1->boundary[i] = MIN(r1->boundary[i], r2->boundary[i]);
+	j = i + t->ndims_alloc;
+	r1->boundary[j] = MAX(r1->boundary[j], r2->boundary[j]);
     }
 
-    return new_rect;
+    for (i = t->ndims; i < t->ndims_alloc; i++) {
+	r1->boundary[i] = 0;
+	j = i + t->ndims_alloc;
+	r1->boundary[j] = 0;
+    }
 }
 
 
@@ -371,7 +400,7 @@ int RTreeCompareRect(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t
     assert(r && s);
 
     for (i = 0; i < t->ndims; i++) {
-	j = i + NUMDIMS;	/* index for high sides */
+	j = i + t->ndims_alloc;	/* index for high sides */
 	if (r->boundary[i] != s->boundary[i] ||
 	    r->boundary[j] != s->boundary[j]) {
 	    return 0;
@@ -382,7 +411,7 @@ int RTreeCompareRect(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t
 
 
 /*-----------------------------------------------------------------------------
-| Decide whether two rectangles overlap.
+| Decide whether two rectangles overlap or touch.
 -----------------------------------------------------------------------------*/
 int RTreeOverlap(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
 {
@@ -391,7 +420,7 @@ int RTreeOverlap(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
     assert(r && s);
 
     for (i = 0; i < t->ndims; i++) {
-	j = i + NUMDIMS;	/* index for high sides */
+	j = i + t->ndims_alloc;	/* index for high sides */
 	if (r->boundary[i] > s->boundary[j] ||
 	    s->boundary[i] > r->boundary[j]) {
 	    return FALSE;
@@ -411,16 +440,16 @@ int RTreeContained(struct RTree_Rect *r, struct RTree_Rect *s, struct RTree *t)
     assert(r && s);
 
     /* undefined rect is contained in any other */
-    if (Undefined(r))
+    if (Undefined(r, t))
 	return TRUE;
 
     /* no rect (except an undefined one) is contained in an undef rect */
-    if (Undefined(s))
+    if (Undefined(s, t))
 	return FALSE;
 
     result = TRUE;
     for (i = 0; i < t->ndims; i++) {
-	j = i + NUMDIMS;	/* index for high sides */
+	j = i + t->ndims_alloc;	/* index for high sides */
 	result = result && r->boundary[i] >= s->boundary[i]
 	    && r->boundary[j] <= s->boundary[j];
     }

lib/vector/rtree/rtree.h

@@ -19,13 +19,11 @@
 #define _R_TREE_H_
 
 #include <stdio.h>
+#include <string.h>
 #include <sys/types.h>
 #include <grass/config.h> /* needed for LFS */
 
 
-#define NUMDIMS	3		/* maximum number of dimensions */
-#define NUMSIDES 2*NUMDIMS
-
 typedef double RectReal;
 
 /*-----------------------------------------------------------------------------
@@ -54,7 +52,7 @@ typedef double RectReal;
 
 struct RTree_Rect
 {
-    RectReal boundary[NUMSIDES];	/* xmin,ymin,...,xmax,ymax,... */
+    RectReal *boundary;	/* xmin,ymin,...,xmax,ymax,... */
 };
 
 struct RTree_Node;               /* node for spatial index */
@@ -72,7 +70,7 @@ struct RTree_Branch              /* branch for spatial index */
     union RTree_Child child;
 };
 
-struct RTree_Node             /* node for spatial index */
+struct RTree_Node       /* node for spatial index */
 {
     int count;          /* number of branches */
     int level;		/* 0 is leaf, others positive */
@@ -96,16 +94,42 @@ typedef int rt_insert_fn(struct RTree_Rect *, union RTree_Child, int, struct RTr
 typedef int rt_delete_fn(struct RTree_Rect *, union RTree_Child, struct RTree *);
 typedef int rt_valid_child_fn(union RTree_Child *);
 
+/* temp vars for each tree */
+/* stacks used for non-recursive insertion/deletion */
+/* stack for file-based index */
+struct fstack
+{
+    struct RTree_Node sn;	/* stack node */
+    int branch_id;		/* branch number to follow down */
+    off_t pos;			/* file position of stack node */
+};
+/* stack for memory-based index */
+struct mstack
+{
+    struct RTree_Node *sn;	/* stack node pointer */
+    int branch_id;		/* branch number to follow down */
+};
+
+/* temp vars for node splitting */
+struct RTree_PartitionVars {
+    int partition[MAXCARD + 1];
+    int total, minfill;
+    int taken[MAXCARD + 1];
+    int count[2];
+    struct RTree_Rect cover[2];
+    RectReal area[2];
+};
+
 struct RTree
 {
     /* RTree setup info */
     int fd;                 	  /* file descriptor */
     unsigned char ndims;    	  /* number of dimensions */
     unsigned char nsides;   	  /* number of sides = 2 * ndims */
-    unsigned char nsides_alloc;   /* number of sides = 2 * ndims */
+    unsigned char ndims_alloc;    	  /* number of dimensions allocated */
+    unsigned char nsides_alloc;   /* number of sides allocated = 2 * ndims allocated */
     int nodesize;           	  /* node size in bytes */
     int branchsize;         	  /* branch size in bytes */
-    int rectsize;           	  /* rectangle size in bytes */
 
     /* RTree info, useful to calculate space requirements */
     int n_nodes;            /* number of nodes */
@@ -149,6 +173,22 @@ struct RTree
     
     struct RTree_Node *root;     /* pointer to root node */
 
+    /* internal variables, specific for each tree,
+     * allocated with tree initialization */
+    /* stacks for tree traversal */
+    struct mstack ms[MAXLEVEL];
+    struct fstack fs[MAXLEVEL];
+    
+    /* variables for splitting / forced reinsertion */
+    struct RTree_PartitionVars p;
+    struct RTree_Branch BranchBuf[MAXCARD + 1];
+
+    struct RTree_Branch tmpb1, tmpb2, c;
+    int BranchCount;
+
+    struct RTree_Rect rect_0, rect_1, upperrect, orect;
+    RectReal *center_n;
+
     off_t rootpos;         /* root node position in file */
 };
 
@@ -167,6 +207,8 @@ void RTreeInitNode(struct RTree *, struct RTree_Node *, int);
 void RTreeFreeNode(struct RTree_Node *);
 void RTreeDestroyNode(struct RTree_Node *, int);
 
+void RTreeNewRect(struct RTree_Rect *, struct RTree *);
+
 /* RTree IO */
 size_t RTreeReadNode(struct RTree_Node *, off_t, struct RTree *);
 size_t RTreeWriteNode(struct RTree_Node *, struct RTree *);

lib/vector/rtree/split.c

@@ -17,9 +17,12 @@
 * 
 ***********************************************************************/
 
+#include <stdlib.h>
 #include <stdio.h>
 #include <assert.h>
 #include <float.h>
+/* remove after debug */
+#include <grass/gis.h>
 #include "index.h"
 #include "card.h"
 #include "split.h"
@@ -28,19 +31,11 @@
 #define DBL_MAX 1.797693E308  /* DBL_MAX approximation */
 #endif
 
-struct RTree_Branch BranchBuf[MAXCARD + 1];
-int BranchCount;
-struct RTree_Rect CoverSplit;
-RectReal CoverSplitArea;
-
-/* variables for finding a partition */
-struct RTree_PartitionVars Partitions[1];
-
 /*----------------------------------------------------------------------
 | Load branch buffer with branches from full node plus the extra branch.
 ----------------------------------------------------------------------*/
 static void RTreeGetBranches(struct RTree_Node *n, struct RTree_Branch *b,
-			     struct RTree *t)
+			     RectReal *CoverSplitArea, struct RTree *t)
 {
     int i, maxkids = 0;
 
@@ -49,7 +44,7 @@ static void RTreeGetBranches(struct RTree_Node *n, struct RTree_Branch *b,
 	/* load the branch buffer */
 	for (i = 0; i < maxkids; i++) {
 	    assert(t->valid_child(&(n->branch[i].child)));	/* n should have every entry full */
-	    BranchBuf[i] = n->branch[i];
+	    RTreeCopyBranch(&(t->BranchBuf[i]), &(n->branch[i]), t);
 	}
     }
     else {
@@ -57,20 +52,20 @@ static void RTreeGetBranches(struct RTree_Node *n, struct RTree_Branch *b,
 	/* load the branch buffer */
 	for (i = 0; i < maxkids; i++) {
 	    assert(n->branch[i].child.id);	/* n should have every entry full */
-	    BranchBuf[i] = n->branch[i];
+	    RTreeCopyBranch(&(t->BranchBuf[i]), &(n->branch[i]), t);
 	}
     }
 
-    BranchBuf[maxkids] = *b;
-    BranchCount = maxkids + 1;
+    RTreeCopyBranch(&(t->BranchBuf[maxkids]), b, t);
+    t->BranchCount = maxkids + 1;
 
     if (METHOD == 0) { /* quadratic split */
 	/* calculate rect containing all in the set */
-	CoverSplit = BranchBuf[0].rect;
+	RTreeCopyRect(&(t->orect), &(t->BranchBuf[0].rect), t);
 	for (i = 1; i < maxkids + 1; i++) {
-	    CoverSplit = RTreeCombineRect(&CoverSplit, &BranchBuf[i].rect, t);
+	    RTreeExpandRect(&(t->orect), &(t->BranchBuf[i].rect), t);
 	}
-	CoverSplitArea = RTreeRectSphericalVolume(&CoverSplit, t);
+	*CoverSplitArea = RTreeRectSphericalVolume(&(t->orect), t);
     }
 
     RTreeInitNode(t, n, NODETYPE(n->level, t->fd));
@@ -89,11 +84,10 @@ static void RTreeClassify(int i, int group, struct RTree_PartitionVars *p,
 
     if (METHOD == 0) {
 	if (p->count[group] == 0)
-	    p->cover[group] = BranchBuf[i].rect;
+	    RTreeCopyRect(&(p->cover[group]), &(t->BranchBuf[i].rect), t);
 	else
-	    p->cover[group] =
-		RTreeCombineRect(&BranchBuf[i].rect, &p->cover[group], t);
-	p->area[group] = RTreeRectSphericalVolume(&p->cover[group], t);
+	    RTreeExpandRect(&(p->cover[group]), &(t->BranchBuf[i].rect), t);
+	p->area[group] = RTreeRectSphericalVolume(&(p->cover[group]), t);
     }
     p->count[group]++;
 }
@@ -109,23 +103,22 @@ static void RTreeClassify(int i, int group, struct RTree_PartitionVars *p,
 | Pick the two that waste the most area if covered by a single
 | rectangle.
 ----------------------------------------------------------------------*/
-static void RTreePickSeeds(struct RTree_PartitionVars *p, struct RTree *t)
+static void RTreePickSeeds(struct RTree_PartitionVars *p, RectReal CoverSplitArea, struct RTree *t)
 {
     int i, j, seed0 = 0, seed1 = 0;
     RectReal worst, waste, area[MAXCARD + 1];
 
     for (i = 0; i < p->total; i++)
-	area[i] = RTreeRectSphericalVolume(&BranchBuf[i].rect, t);
+	area[i] = RTreeRectSphericalVolume(&(t->BranchBuf[i]).rect, t);
 
     worst = -CoverSplitArea - 1;
     for (i = 0; i < p->total - 1; i++) {
 	for (j = i + 1; j < p->total; j++) {
-	    struct RTree_Rect one_rect;
 
-	    one_rect = RTreeCombineRect(&BranchBuf[i].rect,
-					&BranchBuf[j].rect, t);
+	    RTreeCombineRect(&(t->BranchBuf[i].rect), &(t->BranchBuf[j].rect),
+	                     &(t->orect), t);
 	    waste =
-		RTreeRectSphericalVolume(&one_rect, t) - area[i] - area[j];
+		RTreeRectSphericalVolume(&(t->orect), t) - area[i] - area[j];
 	    if (waste > worst) {
 		worst = waste;
 		seed0 = i;
@@ -149,22 +142,25 @@ static void RTreeLoadNodes(struct RTree_Node *n, struct RTree_Node *q,
     for (i = 0; i < p->total; i++) {
 	assert(p->partition[i] == 0 || p->partition[i] == 1);
 	if (p->partition[i] == 0)
-	    RTreeAddBranch(&BranchBuf[i], n, NULL, NULL, NULL, NULL, t);
+	    RTreeAddBranch(&(t->BranchBuf[i]), n, NULL, NULL, NULL, NULL, t);
 	else if (p->partition[i] == 1)
-	    RTreeAddBranch(&BranchBuf[i], q, NULL, NULL, NULL, NULL, t);
+	    RTreeAddBranch(&(t->BranchBuf[i]), q, NULL, NULL, NULL, NULL, t);
     }
 }
 
 /*----------------------------------------------------------------------
 | Initialize a PartitionVars structure.
 ----------------------------------------------------------------------*/
-void RTreeInitPVars(struct RTree_PartitionVars *p, int maxrects, int minfill)
+void RTreeInitPVars(struct RTree_PartitionVars *p, int maxrects, int minfill, struct RTree *t)
 {
     int i;
 
     p->count[0] = p->count[1] = 0;
-    p->cover[0] = p->cover[1] = RTreeNullRect();
-    p->area[0] = p->area[1] = (RectReal) 0;
+    if (METHOD == 0) {
+	RTreeNullRect(&(p->cover[0]), t);
+	RTreeNullRect(&(p->cover[1]), t);
+	p->area[0] = p->area[1] = (RectReal) 0;
+    }
     p->total = maxrects;
     p->minfill = minfill;
     for (i = 0; i < maxrects; i++) {
@@ -177,7 +173,7 @@ void RTreeInitPVars(struct RTree_PartitionVars *p, int maxrects, int minfill)
 | Print out data for a partition from PartitionVars struct.
 | Unused, for debugging only
 ----------------------------------------------------------------------*/
-static void RTreePrintPVars(struct RTree_PartitionVars *p)
+static void RTreePrintPVars(struct RTree_PartitionVars *p, struct RTree *t, RectReal CoverSplitArea)
 {
     int i;
 
@@ -206,10 +202,10 @@ static void RTreePrintPVars(struct RTree_PartitionVars *p)
 		(float)CoverSplitArea / (p->area[0] + p->area[1]));
     }
     fprintf(stdout, "cover[0]:\n");
-    RTreePrintRect(&p->cover[0], 0);
+    RTreePrintRect(&p->cover[0], 0, t);
 
     fprintf(stdout, "cover[1]:\n");
-    RTreePrintRect(&p->cover[1], 0);
+    RTreePrintRect(&p->cover[1], 0, t);
 }
 
 /*----------------------------------------------------------------------
@@ -227,14 +223,18 @@ static void RTreePrintPVars(struct RTree_PartitionVars *p)
 | These last are the ones that can go in either group most easily.
 ----------------------------------------------------------------------*/
 static void RTreeMethodZero(struct RTree_PartitionVars *p, int minfill,
-			    struct RTree *t)
+			    RectReal CoverSplitArea, struct RTree *t)
 {
     int i;
     RectReal biggestDiff;
     int group, chosen = 0, betterGroup = 0;
+    struct RTree_Rect *r, *rect_0, *rect_1;
 
-    RTreeInitPVars(p, BranchCount, minfill);
-    RTreePickSeeds(p, t);
+    RTreeInitPVars(p, t->BranchCount, minfill, t);
+    RTreePickSeeds(p, CoverSplitArea, t);
+    
+    rect_0 = &(t->rect_0);
+    rect_1 = &(t->rect_1);
 
     while (p->count[0] + p->count[1] < p->total
 	   && p->count[0] < p->total - p->minfill
@@ -242,14 +242,13 @@ static void RTreeMethodZero(struct RTree_PartitionVars *p, int minfill,
 	biggestDiff = (RectReal) - 1.;
 	for (i = 0; i < p->total; i++) {
 	    if (!p->taken[i]) {
-		struct RTree_Rect *r, rect_0, rect_1;
 		RectReal growth0, growth1, diff;
 
-		r = &BranchBuf[i].rect;
-		rect_0 = RTreeCombineRect(r, &p->cover[0], t);
-		rect_1 = RTreeCombineRect(r, &p->cover[1], t);
-		growth0 = RTreeRectSphericalVolume(&rect_0, t) - p->area[0];
-		growth1 = RTreeRectSphericalVolume(&rect_1, t) - p->area[1];
+		r = &(t->BranchBuf[i].rect);
+		RTreeCombineRect(r, &(p->cover[0]), rect_0, t);
+		RTreeCombineRect(r, &(p->cover[1]), rect_1, t);
+		growth0 = RTreeRectSphericalVolume(rect_0, t) - p->area[0];
+		growth1 = RTreeRectSphericalVolume(rect_1, t) - p->area[1];
 		diff = growth1 - growth0;
 		if (diff >= 0)
 		    group = 0;
@@ -298,13 +297,11 @@ static void RTreeMethodZero(struct RTree_PartitionVars *p, int minfill,
 /*----------------------------------------------------------------------
 | swap branches
 ----------------------------------------------------------------------*/
-static void RTreeSwapBranches(struct RTree_Branch *a, struct RTree_Branch *b)
+static void RTreeSwapBranches(struct RTree_Branch *a, struct RTree_Branch *b, struct RTree *t)
 {
-    struct RTree_Branch c;
-
-    c = *a;
-    *a = *b;
-    *b = c;
+    RTreeCopyBranch(&(t->c), a, t);
+    RTreeCopyBranch(a, b, t);
+    RTreeCopyBranch(b, &(t->c), t);
 }
 
 /*----------------------------------------------------------------------
@@ -327,12 +324,12 @@ static int RTreeCompareBranches(struct RTree_Branch *a, struct RTree_Branch *b,
 /*----------------------------------------------------------------------
 | check if BranchBuf is sorted along given axis (dimension)
 ----------------------------------------------------------------------*/
-static int RTreeBranchBufIsSorted(int first, int last, int side)
+static int RTreeBranchBufIsSorted(int first, int last, int side, struct RTree *t)
 {
     int i;
 
     for (i = first; i < last; i++) {
-	if (RTreeCompareBranches(&(BranchBuf[i]), &(BranchBuf[i + 1]), side)
+	if (RTreeCompareBranches(&(t->BranchBuf[i]), &(t->BranchBuf[i + 1]), side)
 	    == 1)
 	    return 0;
     }
@@ -343,21 +340,21 @@ static int RTreeBranchBufIsSorted(int first, int last, int side)
 /*----------------------------------------------------------------------
 | partition BranchBuf for quicksort along given axis (dimension)
 ----------------------------------------------------------------------*/
-static int RTreePartitionBranchBuf(int first, int last, int side)
+static int RTreePartitionBranchBuf(int first, int last, int side, struct RTree *t)
 {
     int pivot, mid = (first + last) / 2;
     int larger, smaller;
 
     if (last - first == 1) {	/* only two items in list */
 	if (RTreeCompareBranches
-	    (&(BranchBuf[first]), &(BranchBuf[last]), side) == 1) {
-	    RTreeSwapBranches(&(BranchBuf[first]), &(BranchBuf[last]));
+	    (&(t->BranchBuf[first]), &(t->BranchBuf[last]), side) == 1) {
+	    RTreeSwapBranches(&(t->BranchBuf[first]), &(t->BranchBuf[last]), t);
 	}
 	return last;
     }
 
     /* larger of two */
-    if (RTreeCompareBranches(&(BranchBuf[first]), &(BranchBuf[mid]), side)
+    if (RTreeCompareBranches(&(t->BranchBuf[first]), &(t->BranchBuf[mid]), side)
 	== 1) {
 	larger = pivot = first;
 	smaller = mid;
@@ -367,11 +364,11 @@ static int RTreePartitionBranchBuf(int first, int last, int side)
 	smaller = first;
     }
 
-    if (RTreeCompareBranches(&(BranchBuf[larger]), &(BranchBuf[last]), side)
+    if (RTreeCompareBranches(&(t->BranchBuf[larger]), &(t->BranchBuf[last]), side)
 	== 1) {
 	/* larger is largest, get the larger of smaller and last */
 	if (RTreeCompareBranches
-	    (&(BranchBuf[smaller]), &(BranchBuf[last]), side) == 1) {
+	    (&(t->BranchBuf[smaller]), &(t->BranchBuf[last]), side) == 1) {
 	    pivot = smaller;
 	}
 	else {
@@ -380,16 +377,16 @@ static int RTreePartitionBranchBuf(int first, int last, int side)
     }
 
     if (pivot != last) {
-	RTreeSwapBranches(&(BranchBuf[pivot]), &(BranchBuf[last]));
+	RTreeSwapBranches(&(t->BranchBuf[pivot]), &(t->BranchBuf[last]), t);
     }
 
     pivot = first;
 
     while (first < last) {
 	if (RTreeCompareBranches
-	    (&(BranchBuf[first]), &(BranchBuf[last]), side) != 1) {
+	    (&(t->BranchBuf[first]), &(t->BranchBuf[last]), side) != 1) {
 	    if (pivot != first) {
-		RTreeSwapBranches(&(BranchBuf[pivot]), &(BranchBuf[first]));
+		RTreeSwapBranches(&(t->BranchBuf[pivot]), &(t->BranchBuf[first]), t);
 	    }
 	    pivot++;
 	}
@@ -397,7 +394,7 @@ static int RTreePartitionBranchBuf(int first, int last, int side)
     }
 
     if (pivot != last) {
-	RTreeSwapBranches(&(BranchBuf[pivot]), &(BranchBuf[last]));
+	RTreeSwapBranches(&(t->BranchBuf[pivot]), &(t->BranchBuf[last]), t);
     }
 
     return pivot;
@@ -406,13 +403,13 @@ static int RTreePartitionBranchBuf(int first, int last, int side)
 /*----------------------------------------------------------------------
 | quicksort BranchBuf along given side
 ----------------------------------------------------------------------*/
-static void RTreeQuicksortBranchBuf(int side)
+static void RTreeQuicksortBranchBuf(int side, struct RTree *t)
 {
     int pivot, first, last;
     int s_first[MAXCARD + 1], s_last[MAXCARD + 1], stacksize;
 
     s_first[0] = 0;
-    s_last[0] = BranchCount - 1;
+    s_last[0] = t->BranchCount - 1;
 
     stacksize = 1;
 
@@ -422,9 +419,9 @@ static void RTreeQuicksortBranchBuf(int side)
 	first = s_first[stacksize];
 	last = s_last[stacksize];
 	if (first < last) {
-	    if (!RTreeBranchBufIsSorted(first, last, side)) {
+	    if (!RTreeBranchBufIsSorted(first, last, side, t)) {
 
-		pivot = RTreePartitionBranchBuf(first, last, side);
+		pivot = RTreePartitionBranchBuf(first, last, side, t);
 
 		s_first[stacksize] = first;
 		s_last[stacksize] = pivot - 1;
@@ -454,14 +451,29 @@ static void RTreeMethodOne(struct RTree_PartitionVars *p, int minfill,
                            int maxkids, struct RTree *t)
 {
     int i, j, k, l, s;
-    int axis = 0, best_axis = 0, side = 0, best_side[NUMDIMS];
-    int best_cut[NUMDIMS];
+    int axis = 0, best_axis = 0, side = 0;
     RectReal margin, smallest_margin = 0;
-    struct RTree_Rect *r1, *r2, testrect1, testrect2, upperrect, orect;
+    struct RTree_Rect *r1, *r2;
+    struct RTree_Rect *rect_0, *rect_1, *orect, *upperrect;
     int minfill1 = minfill - 1;
     RectReal overlap, vol, smallest_overlap = -1, smallest_vol = -1;
 
-    RTreeInitPVars(p, BranchCount, minfill);
+    static int *best_cut = NULL, *best_side = NULL;
+    static int one_init = 0;
+    
+    if (!one_init) {
+	best_cut = (int *)malloc(t->ndims_alloc * sizeof(int));
+	best_side = (int *)malloc(t->ndims_alloc * sizeof(int));
+	one_init = 1;
+    }
+
+    rect_0 = &(t->rect_0);
+    rect_1 = &(t->rect_1);
+    orect = &(t->orect);
+    upperrect = &(t->upperrect);
+
+    RTreeInitPVars(p, t->BranchCount, minfill, t);
+    RTreeInitRect(orect, t);
 
     margin = DBL_MAX;
 
@@ -478,38 +490,38 @@ static void RTreeMethodOne(struct RTree_PartitionVars *p, int minfill,
 
 	/* first lower then upper bounds for each axis */
 	for (s = 0; s < 2; s++) {
-	    RTreeQuicksortBranchBuf(i + s * NUMDIMS);
+	    RTreeQuicksortBranchBuf(i + s * t->ndims_alloc, t);
 
 	    side = s;
 
-	    testrect1 = BranchBuf[0].rect;
-	    upperrect = BranchBuf[maxkids].rect;
+	    RTreeCopyRect(rect_0, &(t->BranchBuf[0].rect), t);
+	    RTreeCopyRect(upperrect, &(t->BranchBuf[maxkids].rect), t);
 
 	    for (j = 1; j < minfill1; j++) {
-		r1 = &BranchBuf[j].rect;
-		testrect1 = RTreeCombineRect(&testrect1, r1, t);
-		r2 = &BranchBuf[maxkids - j].rect;
-		upperrect = RTreeCombineRect(&upperrect, r2, t);
+		r1 = &(t->BranchBuf[j].rect);
+		RTreeExpandRect(rect_0, r1, t);
+		r2 = &(t->BranchBuf[maxkids - j].rect);
+		RTreeExpandRect(upperrect, r2, t);
 	    }
-	    r2 = &BranchBuf[maxkids - minfill1].rect;
-	    upperrect = RTreeCombineRect(&upperrect, r2, t);
+	    r2 = &(t->BranchBuf[maxkids - minfill1].rect);
+	    RTreeExpandRect(upperrect, r2, t);
 
 	    /* check distributions for this axis, adhere the minimum node fill */
-	    for (j = minfill1; j < BranchCount - minfill; j++) {
+	    for (j = minfill1; j < t->BranchCount - minfill; j++) {
 
-		r1 = &BranchBuf[j].rect;
-		testrect1 = RTreeCombineRect(&testrect1, r1, t);
+		r1 = &(t->BranchBuf[j].rect);
+		RTreeExpandRect(rect_0, r1, t);
 
-		testrect2 = upperrect;
-		for (k = j + 1; k < BranchCount - minfill; k++) {
-		    r2 = &BranchBuf[k].rect;
-		    testrect2 = RTreeCombineRect(&testrect2, r2, t);
+		RTreeCopyRect(rect_1, upperrect, t);
+		for (k = j + 1; k < t->BranchCount - minfill; k++) {
+		    r2 = &(t->BranchBuf[k].rect);
+		    RTreeExpandRect(rect_1, r2, t);
 		}
 
 		/* the margin is the sum of the lengths of the edges of a rectangle */
 		margin =
-		    RTreeRectMargin(&testrect1, t) +
-		    RTreeRectMargin(&testrect2, t);
+		    RTreeRectMargin(rect_0, t) +
+		    RTreeRectMargin(rect_1, t);
 
 		/* remember best axis */
 		if (margin <= smallest_margin) {
@@ -524,31 +536,31 @@ static void RTreeMethodOne(struct RTree_PartitionVars *p, int minfill,
 
 		for (k = 0; k < t->ndims; k++) {
 		    /* no overlap */
-		    if (testrect1.boundary[k] > testrect2.boundary[k + NUMDIMS] ||
-			testrect1.boundary[k + NUMDIMS] < testrect2.boundary[k]) {
+		    if (rect_0->boundary[k] > rect_1->boundary[k + t->ndims_alloc] ||
+			rect_0->boundary[k + t->ndims_alloc] < rect_1->boundary[k]) {
 			overlap = 0;
 			break;
 		    }
 		    /* get overlap */
 		    else {
-			if (testrect1.boundary[k] > testrect2.boundary[k])
-			    orect.boundary[k] = testrect1.boundary[k];
+			if (rect_0->boundary[k] > rect_1->boundary[k])
+			    orect->boundary[k] = rect_0->boundary[k];
 			else
-			    orect.boundary[k] = testrect2.boundary[k];
+			    orect->boundary[k] = rect_1->boundary[k];
 
-			l = k + NUMDIMS;
-			if (testrect1.boundary[l] < testrect2.boundary[l])
-			    orect.boundary[l] = testrect1.boundary[l];
+			l = k + t->ndims_alloc;
+			if (rect_0->boundary[l] < rect_1->boundary[l])
+			    orect->boundary[l] = rect_0->boundary[l];
 			else
-			    orect.boundary[l] = testrect2.boundary[l];
+			    orect->boundary[l] = rect_1->boundary[l];
 		    }
 		}
 		if (overlap)
-		    overlap = RTreeRectVolume(&orect, t);
+		    overlap = RTreeRectVolume(orect, t);
 
 		vol =
-		    RTreeRectVolume(&testrect1, t) +
-		    RTreeRectVolume(&testrect2, t);
+		    RTreeRectVolume(rect_0, t) +
+		    RTreeRectVolume(rect_1, t);
 
 		/* get best cut for this axis */
 		if (overlap <= smallest_overlap) {
@@ -571,14 +583,14 @@ static void RTreeMethodOne(struct RTree_PartitionVars *p, int minfill,
 
     /* Use best distribution to classify branches */
     if (best_axis != axis || best_side[best_axis] != side)
-	RTreeQuicksortBranchBuf(best_axis + best_side[best_axis] * NUMDIMS);
+	RTreeQuicksortBranchBuf(best_axis + best_side[best_axis] * t->ndims_alloc, t);
 
     best_cut[best_axis]++;
 
     for (i = 0; i < best_cut[best_axis]; i++)
 	RTreeClassify(i, 0, p, t);
 
-    for (i = best_cut[best_axis]; i < BranchCount; i++)
+    for (i = best_cut[best_axis]; i < t->BranchCount; i++)
 	RTreeClassify(i, 1, p, t);
 
     assert(p->count[0] + p->count[1] == p->total);
@@ -595,19 +607,20 @@ void RTreeSplitNode(struct RTree_Node *n, struct RTree_Branch *b,
                     struct RTree_Node *nn, struct RTree *t)
 {
     struct RTree_PartitionVars *p;
+    RectReal CoverSplitArea;
     int level;
 
     /* load all the branches into a buffer, initialize old node */
     level = n->level;
-    RTreeGetBranches(n, b, t);
+    RTreeGetBranches(n, b, &CoverSplitArea, t);
 
     /* find partition */
-    p = &Partitions[0];
+    p = &(t->p);
     
     if (METHOD == 1) /* R* split */
-	RTreeMethodOne(p, MINFILL(level, t), MAXKIDS(level, t), t);
+	RTreeMethodOne(&(t->p), MINFILL(level, t), MAXKIDS(level, t), t);
     else
-	RTreeMethodZero(p, MINFILL(level, t), t);
+	RTreeMethodZero(&(t->p), MINFILL(level, t), CoverSplitArea, t);
 
     /*
      * put branches from buffer into 2 nodes

lib/vector/rtree/split.h

@@ -24,17 +24,4 @@
 /* METHOD 1: R*-Tree split */
 #define METHOD 1
 
-struct RTree_PartitionVars {
-    int partition[MAXCARD + 1];
-    int total, minfill;
-    int taken[MAXCARD + 1];
-    int count[2];
-    struct RTree_Rect cover[2];
-    RectReal area[2];
-};
-
-
-/* variables for finding a partition */
-extern struct RTree_PartitionVars Partitions[1];
-
-extern void RTreeInitPVars(struct RTree_PartitionVars *, int, int);
+void RTreeInitPVars(struct RTree_PartitionVars *, int, int, struct RTree *);