new rtree

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

lib/vector/rtree/Makefile

@@ -12,7 +12,7 @@ ifneq ($(USE_LARGEFILES),)
 endif
 
 RTLINC = $(ARCH_INCDIR)/rtree
-HEADERS := $(RTLINC)/card.h $(RTLINC)/index.h $(RTLINC)/split_q.h \
+HEADERS := $(RTLINC)/card.h $(RTLINC)/index.h $(RTLINC)/split.h \
 	$(ARCH_INCDIR)/rtree.h
 
 default: headers

lib/vector/rtree/Makefile.alone

@@ -8,7 +8,7 @@ OBJECTS = card.o \
 	node.o \
 	rect.o \
 	sphvol.o \
-	split_q.o 
+	split.o 
 
 all: librtree.a test
 

lib/vector/rtree/card.c

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -29,19 +30,19 @@ static int set_max(int *which, int new_max)
     return 1;
 }
 
-int RTreeSetNodeMax(int new_max)
+int RTreeSetNodeMax(int new_max, struct RTree *t)
 {
-    return set_max(&NODECARD, new_max);
+    return set_max(&(t->nodecard), new_max);
 }
-int RTreeSetLeafMax(int new_max)
+int RTreeSetLeafMax(int new_max, struct RTree *t)
 {
-    return set_max(&LEAFCARD, new_max);
+    return set_max(&(t->leafcard), new_max);
 }
-int RTreeGetNodeMax(void)
+int RTreeGetNodeMax(struct RTree *t)
 {
-    return NODECARD;
+    return t->nodecard;
 }
-int RTreeGetLeafMax(void)
+int RTreeGetLeafMax(struct RTree *t)
 {
-    return LEAFCARD;
+    return t->leafcard;
 }

lib/vector/rtree/card.h

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS

lib/vector/rtree/docs/README.grass

@@ -1,9 +1,10 @@
-This library was taken from:
+This library was originally taken from:
 http://www.superliminal.com/sources/RTree.zip
 described on:
 http://www.superliminal.com/sources/sources.htm
 Copy of this file is in sources.htm
-
+The current R*-tree implementation is based on
+http://dbs.mathematik.uni-marburg.de/publications/myPapers/1990/BKSS90.pdf
 
 Changes done in GRASS:
 - files converted to unix line ends
@@ -20,4 +21,4 @@ Changes done in GRASS:
 - type RectReal: float -> double 
 - log_pi set to constant
 - added GRASS headers to all *.[hc] files
-
+- changed to R*-tree, near complete rewrite

lib/vector/rtree/docs/README.txt

@@ -1,8 +1,10 @@
 
-This directory contains a C implementation of the R-tree data structure.
-The implementation is originally from the test code of the authors, and
-later ported to ANSI C on a variety of platforms by Daniel Green
-(dgreen@superliminal.com).
+This directory contains a C implementation of the R*-tree data structure.
+The implementation is originally from the RTree test code of Toni 
+Guttmann, and later ported to ANSI C on a variety of platforms by 
+Daniel Green (dgreen@superliminal.com). 
+Later on it was converted to an R*-tree by Markus Metz based on the article
+http://dbs.mathematik.uni-marburg.de/publications/myPapers/1990/BKSS90.pdf.
 
 Paul Brooke (pbrooke@mindscape.com) discovered an interesting anomaly in
 the original algorithm which uses the rectangular volumes of nodes as the

lib/vector/rtree/docs/test.c

@@ -40,7 +40,7 @@ int MySearchCallback(int id, void *arg)
 
 int main()
 {
-    struct Node *root = RTreeNewIndex();
+    struct RTree *rtree = RTreeNewIndex(2);
     int i, nhits;
 
     fprintf(stdout, "nrects = %d\n", nrects);
@@ -55,8 +55,8 @@ int main()
      * parameter 4 is always zero which means to add from the root.
      */
     for (i = 0; i < nrects; i++)
-	RTreeInsertRect(&rects[i], i + 1, &root, 0);	/* i+1 is rect ID. Note: root can change */
-    nhits = RTreeSearch(root, &search_rect, MySearchCallback, 0);
+	RTreeInsertRect(&rects[i], i + 1, rtree);	/* i+1 is rect ID. */
+    nhits = RTreeSearch(rtree, &search_rect, MySearchCallback, 0);
     fprintf(stdout, "Search resulted in %d hits\n", nhits);
 
     return 0;

lib/vector/rtree/gammavol.c

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -16,7 +17,6 @@
 *****************************************************************************/
 #include <stdio.h>
 #include <math.h>
-#include <grass/gis.h>
 
 #ifndef ABS
 #	define ABS(a) ((a) > 0 ? (a) : -(a))

lib/vector/rtree/index.c

@@ -5,178 +5,331 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *****************************************************************************/
 
-#include <stdio.h>
 #include <stdlib.h>
-#include "assert.h"
+#include <assert.h>
 #include "index.h"
 #include "card.h"
 
-/* Make a new index, empty.  Consists of a single node. */
-struct Node *RTreeNewIndex(void)
+/* stack used for non-recursive insertion/deletion */
+struct stack
 {
-    struct Node *x;
+    struct Node *sn;		/* stack node */
+    int branch_id;		/* branch no to follow down */
+};
 
-    x = RTreeNewNode();
-    x->level = 0;		/* leaf */
-    return x;
+/* 
+ * Make a new index, empty.
+ * ndims number of dimensions
+ * returns pointer to RTree structure
+ */
+struct RTree *RTreeNewIndex(int ndims)
+{
+    struct RTree *new_rtree;
+    struct Node *n;
+
+    new_rtree = (struct RTree *)malloc(sizeof(struct RTree));
+
+    new_rtree->ndims = ndims;
+    new_rtree->nsides = 2 * ndims;
+
+    new_rtree->nodesize = sizeof(struct Node);
+    new_rtree->branchsize = sizeof(struct Branch);
+    new_rtree->rectsize = sizeof(struct Rect);
+
+    /* nodecard and leafcard can be adjusted, must NOT be larger than MAXCARD */
+    new_rtree->nodecard = MAXCARD;
+    new_rtree->leafcard = MAXCARD;
+    /* NOTE: min fill can be changed if needed. */
+    new_rtree->min_node_fill = (new_rtree->nodecard - 1) / 2;
+    new_rtree->min_leaf_fill = (new_rtree->leafcard - 1) / 2;
+
+    n = RTreeNewNode(new_rtree);
+    new_rtree->n_levels = n->level = 0;	/* leaf */
+    new_rtree->root = n;
+
+    new_rtree->n_nodes = 1;
+    new_rtree->n_leafs = 0;
+
+    return new_rtree;
+}
+
+void RTreeFreeIndex(struct RTree *t)
+{
+    assert(t);
+
+    if (t->root)
+	RTreeDestroyNode(t->root, t->nodecard);
+
+    free(t);
 }
 
 /*
- * Search in an index tree or subtree for all data retangles that
+ * Search in an index tree for all data retangles that
  * overlap the argument rectangle.
  * Return the number of qualifying data rects.
  */
-int RTreeSearch(struct Node *N, struct Rect *R, SearchHitCallback shcb,
+int RTreeSearch(struct RTree *t, struct Rect *r, SearchHitCallback shcb,
 		void *cbarg)
 {
-    register struct Node *n = N;
-    register struct Rect *r = R;	/* NOTE: Suspected bug was R sent in as Node* and cast to Rect* here. */
-
-    /* Fix not yet tested. */
-    register int hitCount = 0;
-    register int i;
+    struct Node *n;
+    int hitCount = 0, found;
+    int i;
+    struct stack s[MAXLEVEL];
+    int top = 0;
 
-    assert(n);
-    assert(n->level >= 0);
     assert(r);
-
-    if (n->level > 0) {		/* this is an internal node in the tree */
-	for (i = 0; i < NODECARD; i++)
-	    if (n->branch[i].child && RTreeOverlap(r, &n->branch[i].rect)) {
-		hitCount += RTreeSearch(n->branch[i].child, r, shcb, cbarg);
+    assert(t);
+
+    /* stack size of t->n_levels + 1 is enough because of depth first search */
+    /* only one node per level on stack at any given time */
+
+    /* add root node position to stack */
+    s[top].sn = t->root;
+    s[top].branch_id = i = 0;
+    n = s[top].sn;
+
+    while (top >= 0) {
+	n = s[top].sn;
+	if (s[top].sn->level > 0) {	/* this is an internal node in the tree */
+	    found = 1;
+	    for (i = s[top].branch_id; i < t->nodecard; i++) {
+		if (s[top].sn->branch[i].child &&
+		    RTreeOverlap(r, &(s[top].sn->branch[i].rect), t)) {
+		    s[top++].branch_id = i + 1;
+		    /* add next node to stack */
+		    s[top].sn = n->branch[i].child;
+		    s[top].branch_id = 0;
+		    found = 0;
+		    break;
+		}
 	    }
-    }
-    else {			/* this is a leaf node */
-
-	for (i = 0; i < LEAFCARD; i++)
-	    if (n->branch[i].child && RTreeOverlap(r, &n->branch[i].rect)) {
-		hitCount++;
-		if (shcb)	/* call the user-provided callback */
-		    if (!shcb((int)n->branch[i].child, cbarg))
-			return hitCount;	/* callback wants to terminate search early */
+	    if (found) {
+		/* nothing else found, go back up */
+		s[top].branch_id = t->nodecard;
+		top--;
+	    }
+	}
+	else {			/* this is a leaf node */
+	    for (i = 0; i < t->leafcard; i++) {
+		if (s[top].sn->branch[i].child &&
+		    RTreeOverlap(r, &(s[top].sn->branch[i].rect), t)) {
+		    hitCount++;
+		    if (shcb) {	/* call the user-provided callback */
+			if (!shcb((int)s[top].sn->branch[i].child, cbarg)) {
+			    /* callback wants to terminate search early */
+			    return hitCount;
+			}
+		    }
+		}
 	    }
+	    top--;
+	}
     }
+
     return hitCount;
 }
 
+/* 
+ * Free ListBranch
+ */
+static void RTreeFreeListBranch(struct ListBranch *p)
+{
+    free(p);
+}
+
 /*
  * Inserts a new data rectangle into the index structure.
- * Recursively descends tree, propagates splits back up.
- * Returns 0 if node was not split.  Old node updated.
- * If node was split, returns 1 and sets the pointer pointed to by
- * new_node to point to the new node.  Old node updated to become one of two.
+ * Non-recursively descends tree.
+ * Returns 0 if node was not split and nothing was removed.
+ * Returns 1 if root node was split. Old node updated to become one of two.
+ * Returns 2 if branches need to be reinserted.  
  * The level argument specifies the number of steps up from the leaf
  * level to insert; e.g. a data rectangle goes in at level = 0.
  */
-static int RTreeInsertRect2(struct Rect *r,
-			    struct Node *child, struct Node *n, struct Node **new_node,
-			    int level)
+static int RTreeInsertRect2(struct Rect *r, struct Node *child, int level,
+			    struct Node **newnode, struct RTree *t,
+			    struct ListBranch **ee, int *overflow)
 {
-    /*
-       register struct Rect *r = R;
-       register int tid = Tid;
-       register struct Node *n = N, **new_node = New_node;
-       register int level = Level;
-     */
-
-    register int i;
+    int i;
     struct Branch b;
-    struct Node *n2;
-
-    assert(r && n && new_node);
-    assert(level >= 0 && level <= n->level);
-
-    /* Still above level for insertion, go down tree recursively */
-    if (n->level > level) {
-	i = RTreePickBranch(r, n);
-	if (!RTreeInsertRect2(r, child, n->branch[i].child, &n2, level)) {
-	    /* child was not split */
-	    n->branch[i].rect = RTreeCombineRect(r, &(n->branch[i].rect));
-	    return 0;
-	}
-	else {			/* child was split */
+    struct Node *n, *n2;
+    struct Rect *cover;
+    struct stack s[MAXLEVEL];
+    int top = 0, down = 0;
+    int result;
 
-	    n->branch[i].rect = RTreeNodeCover(n->branch[i].child);
-	    b.child = n2;
-	    b.rect = RTreeNodeCover(n2);
-	    return RTreeAddBranch(&b, n, new_node);
-	}
+    assert(r && newnode && t);
+
+    /* add root node to stack */
+    s[top].sn = t->root;
+
+    /* go down to level of insertion */
+    while (s[top].sn->level > level) {
+	n = s[top].sn;
+	i = RTreePickBranch(r, n, t);
+	s[top++].branch_id = i;
+	/* add next node to stack */
+	s[top].sn = n->branch[i].child;
     }
 
-    /* Have reached level for insertion. Add rect, split if necessary */
-    else if (n->level == level) {
+    /* Have reached level for insertion. Remove p rectangles or split */
+    if (s[top].sn->level == level) {
 	b.rect = *r;
-	b.child = child;
 	/* child field of leaves contains tid of data record */
-	return RTreeAddBranch(&b, n, new_node);
+	b.child = child;
+	/* add branch, may split node or remove branches */
+	cover = &(s[top - 1].sn->branch[s[top - 1].branch_id].rect);
+	result = RTreeAddBranch(&b, s[top].sn, &n2, ee, cover, overflow, t);
+	/* update node count */
+	if (result == 1) {
+	    t->n_nodes++;
+	}
     }
     else {
 	/* Not supposed to happen */
 	assert(FALSE);
 	return 0;
     }
+
+    /* go back up */
+    while (top) {
+	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);
+	}
+	else if (result == 2) {	/* branches were removed */
+	    /* get node cover of previous node */
+	    s[top].sn->branch[i].rect = RTreeNodeCover(s[down].sn, 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);
+	    /* add new branch for new node previously added by RTreeAddBranch() */
+	    b.child = n2;
+	    b.rect = RTreeNodeCover(b.child, t);
+
+	    /* add branch, may split node or remove branches */
+	    cover = &(s[top - 1].sn->branch[s[top - 1].branch_id].rect);
+	    result =
+		RTreeAddBranch(&b, s[top].sn, &n2, ee, cover, overflow, t);
+
+	    /* update node count */
+	    if (result == 1) {
+		t->n_nodes++;
+	    }
+	}
+    }
+
+    *newnode = n2;
+
+    return result;
 }
 
 /* 
  * Insert a data rectangle into an index structure.
- * RTreeInsertRect provides for splitting the root;
+ * RTreeInsertRect1 provides for splitting the root;
  * returns 1 if root was split, 0 if it was not.
  * The level argument specifies the number of steps up from the leaf
  * level to insert; e.g. a data rectangle goes in at level = 0.
- * RTreeInsertRect2 does the recursion.
+ * RTreeInsertRect2 does the actual insertion.
  */
-int RTreeInsertRect(struct Rect *R, int Tid, struct Node **Root, int Level)
+static int RTreeInsertRect1(struct Rect *r, struct Node *child, int level,
+			    struct RTree *t)
 {
-    assert(Level == 0);
-    return RTreeInsertRect1(R, (struct Node *) Tid, Root, Level);
-}
-
-int RTreeInsertRect1(struct Rect *R, struct Node *Child, struct Node **Root, int Level)
-{
-    register struct Rect *r = R;
-    register struct Node *child = Child;
-    register struct Node **root = Root;
-    register int level = Level;
-    register int i;
-    register struct Node *newroot;
     struct Node *newnode;
+    struct Node *newroot;
     struct Branch b;
+    struct ListBranch *reInsertList = NULL;
+    struct ListBranch *e;
     int result;
-
-    assert(r && root);
-    assert(level >= 0 && level <= (*root)->level);
-    for (i = 0; i < NUMDIMS; i++) {
-	assert(r->boundary[i] <= r->boundary[NUMDIMS + i]);
-    }
-
-    if (RTreeInsertRect2(r, child, *root, &newnode, level)) {	/* root split */
-	newroot = RTreeNewNode();	/* grow a new root, & tree taller */
-	newroot->level = (*root)->level + 1;
-	b.rect = RTreeNodeCover(*root);
-	b.child = *root;
-	RTreeAddBranch(&b, newroot, NULL);
-	b.rect = RTreeNodeCover(newnode);
+    int i, overflow[50];
+
+    /* R*-tree forced reinsertion: for each level only once */
+    for (i = 0; i < 50; i++)
+	overflow[i] = 1;
+
+    result =
+	RTreeInsertRect2(r, child, level, &newnode, t, &reInsertList,
+			 overflow);
+    if (result == 1) {		/* root split */
+	/* grow a new root, & tree taller */
+	newroot = RTreeNewNode(t);
+	newroot->level = ++t->n_levels;
+	/* branch for old root */
+	b.rect = RTreeNodeCover(t->root, t);
+	b.child = t->root;
+	RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+	/* branch for new node created by RTreeInsertRect2() */
+	b.rect = RTreeNodeCover(newnode, t);
 	b.child = newnode;
-	RTreeAddBranch(&b, newroot, NULL);
-	*root = newroot;
-	result = 1;
+	RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+	/* set new root node */
+	t->root = newroot;
+	t->n_nodes++;
+    }
+    else if (result == 2) {	/* branches were removed */
+	while (reInsertList) {
+	    /* get next branch in list */
+	    b = reInsertList->b;
+	    level = reInsertList->level;
+	    e = reInsertList;
+	    reInsertList = reInsertList->next;
+	    RTreeFreeListBranch(e);
+	    /* reinsert branches */
+	    result =
+		RTreeInsertRect2(&(b.rect), b.child, level, &newnode, t,
+				 &reInsertList, overflow);
+
+	    if (result == 1) {	/* root split */
+		/* grow a new root, & tree taller */
+		newroot = RTreeNewNode(t);
+		newroot->level = ++t->n_levels;
+		/* branch for old root */
+		b.rect = RTreeNodeCover(t->root, t);
+		b.child = t->root;
+		RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+		/* branch for new node created by RTreeInsertRect2() */
+		b.rect = RTreeNodeCover(newnode, t);
+		b.child = newnode;
+		RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+		/* set new root node */
+		t->root = newroot;
+		t->n_nodes++;
+	    }
+	}
     }
-    else
-	result = 0;
 
     return result;
 }
 
+/* 
+ * Insert a data rectangle into an RTree index structure.
+ * r pointer to rectangle
+ * tid data id stored with rectangle, must be > 0
+ * t RTree where rectangle should be inserted
+ */
+int RTreeInsertRect(struct Rect *r, int tid, struct RTree *t)
+{
+    assert(r && t);
+
+    t->n_leafs++;
+
+    return RTreeInsertRect1(r, (struct Node *)tid, 0, t);
+}
+
 /*
  * Allocate space for a node in the list used in DeletRect to
  * store Nodes that are too empty.
@@ -190,14 +343,13 @@ static struct ListNode *RTreeNewListNode(void)
 static void RTreeFreeListNode(struct ListNode *p)
 {
     free(p);
-    /* delete(p); */
 }
 
 /* 
  * Add a node to the reinsertion list.  All its branches will later
  * be reinserted into the index structure.
  */
-static void RTreeReInsert(struct Node *n, struct ListNode **ee)
+static void RTreeReInsertNode(struct Node *n, struct ListNode **ee)
 {
     register struct ListNode *l;
 
@@ -209,94 +361,123 @@ static void RTreeReInsert(struct Node *n, struct ListNode **ee)
 
 /*
  * Delete a rectangle from non-root part of an index structure.
- * Called by RTreeDeleteRect.  Descends tree recursively,
+ * Called by RTreeDeleteRect.  Descends tree non-recursively,
  * merges branches on the way back up.
  * Returns 1 if record not found, 0 if success.
  */
 static int
-RTreeDeleteRect2(struct Rect *R, struct Node *Child, struct Node *N,
-		 struct ListNode **Ee)
+RTreeDeleteRect2(struct Rect *r, struct Node *child, struct RTree *t,
+		 struct ListNode **ee)
 {
-    register struct Rect *r = R;
-    register struct Node *child = Child;
-    register struct Node *n = N;
-    register struct ListNode **ee = Ee;
-    register int i;
-
-    assert(r && n && ee);
-    assert(child);
-    assert(n->level >= 0);
-
-    if (n->level > 0) {		/* not a leaf node */
-	for (i = 0; i < NODECARD; i++) {
-	    if (n->branch[i].child && RTreeOverlap(r, &(n->branch[i].rect))) {
-		if (!RTreeDeleteRect2(r, child, n->branch[i].child, ee)) {
-		    if (n->branch[i].child->count >= MinNodeFill) {
-			n->branch[i].rect =
-			    RTreeNodeCover(n->branch[i].child);
-		    }
-		    else {
-			/* not enough entries in child, eliminate child node */
-			RTreeReInsert(n->branch[i].child, ee);
-			RTreeDisconnectBranch(n, i);
-		    }
-		    return 0;
+    int i, notfound = 1;
+    struct Node *n;
+    struct stack s[MAXLEVEL];
+    int top = 0, down = 0;
+    int minfill;
+
+    assert(r && ee && t);
+
+    /* add root node position to stack */
+    s[top].sn = t->root;
+    s[top].branch_id = 0;
+    n = s[top].sn;
+
+    while (notfound) {
+	/* go down to level 0, remember path */
+	if (s[top].sn->level > 0) {
+	    n = s[top].sn;
+	    for (i = s[top].branch_id; i < t->nodecard; i++) {
+		if (n->branch[i].child &&
+		    RTreeOverlap(r, &(n->branch[i].rect), t)) {
+		    s[top++].branch_id = i + 1;
+		    /* add next node to stack */
+		    s[top].sn = n->branch[i].child;
+		    s[top].branch_id = 0;
+
+		    notfound = 0;
+		    break;
 		}
 	    }
+	    if (notfound) {
+		/* nothing else found, go back up */
+		s[top].branch_id = t->nodecard;
+		top--;
+	    }
+	    else		/* found a way down but not yet the item */
+		notfound = 1;
+	}
+	else {
+	    for (i = 0; i < t->leafcard; i++) {
+		if (s[top].sn->branch[i].child && s[top].sn->branch[i].child == child) {	/* found item */
+		    RTreeDisconnectBranch(s[top].sn, i, t);
+		    t->n_leafs--;
+		    notfound = 0;
+		    break;
+		}
+	    }
+	    if (notfound)	/* continue searching */
+		top--;
 	}
-	return 1;
     }
-    else {			/* a leaf node */
 
-	for (i = 0; i < LEAFCARD; i++) {
-	    if (n->branch[i].child &&
-		n->branch[i].child == child) {
-		RTreeDisconnectBranch(n, i);
-		return 0;
-	    }
+    if (notfound) {
+	return notfound;
+    }
+
+    /* go back up */
+    while (top) {
+	down = top;
+	top--;
+	n = s[top].sn;
+	i = s[top].branch_id - 1;
+	assert(s[down].sn->level == s[top].sn->level - 1);
+
+	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);
+	}
+	else {
+	    /* not enough entries in child, eliminate child node */
+	    RTreeReInsertNode(s[top].sn->branch[i].child, ee);
+	    RTreeDisconnectBranch(s[top].sn, i, t);
 	}
-	return 1;
     }
+
+    return notfound;
 }
 
 /*
+ * should be called by RTreeDeleteRect() only
+ * 
  * Delete a data rectangle from an index structure.
- * Pass in a pointer to a Rect, the tid of the record, ptr to ptr to root node.
+ * Pass in a pointer to a Rect, the tid of the record, ptr RTree.
  * Returns 1 if record not found, 0 if success.
- * RTreeDeleteRect provides for eliminating the root.
+ * RTreeDeleteRect1 provides for eliminating the root.
  */
-int RTreeDeleteRect(struct Rect *R, int Tid, struct Node **Nn)
+static int RTreeDeleteRect1(struct Rect *r, struct Node *child,
+			    struct RTree *t)
 {
-    /* wrapper not really needed, but restricts compile warnings to rtree lib */
-    /* this way it's easier to fix if necessary? */
-    return RTreeDeleteRect1(R, (struct Node *) Tid, Nn);
-}
-
-int RTreeDeleteRect1(struct Rect *R, struct Node *Child, struct Node **Nn)
-{
-    register struct Rect *r = R;
-    register struct Node *child = Child;
-    register struct Node **nn = Nn;
-    register int i;
-    struct Node *tmp_nptr = NULL;
+    int i, maxkids;
+    struct Node *n;
     struct ListNode *reInsertList = NULL;
-    register struct ListNode *e;
+    struct ListNode *e;
 
-    assert(r && nn);
-    assert(*nn);
-    assert(child);
+    assert(r);
+    assert(t);
 
-    if (!RTreeDeleteRect2(r, child, *nn, &reInsertList)) {
+    if (!RTreeDeleteRect2(r, child, t, &reInsertList)) {
 	/* found and deleted a data item */
 
 	/* reinsert any branches from eliminated nodes */
 	while (reInsertList) {
-	    tmp_nptr = reInsertList->node;
-	    for (i = 0; i < MAXKIDS(tmp_nptr); i++) {
-		if (tmp_nptr->branch[i].child) {
-		    RTreeInsertRect1(&(tmp_nptr->branch[i].rect),
-				    tmp_nptr->branch[i].child,
-				    nn, tmp_nptr->level);
+	    t->n_nodes--;
+	    n = reInsertList->node;
+	    maxkids = (n->level > 0 ? t->nodecard : t->leafcard);
+	    for (i = 0; i < maxkids; i++) {
+		if (n->branch[i].child) {
+		    RTreeInsertRect1(&(n->branch[i].rect),
+				     n->branch[i].child, n->level, t);
 		}
 	    }
 	    e = reInsertList;
@@ -306,15 +487,15 @@ int RTreeDeleteRect1(struct Rect *R, struct Node *Child, struct Node **Nn)
 	}
 
 	/* check for redundant root (not leaf, 1 child) and eliminate */
-	if ((*nn)->count == 1 && (*nn)->level > 0) {
-	    for (i = 0; i < NODECARD; i++) {
-		tmp_nptr = (*nn)->branch[i].child;
-		if (tmp_nptr)
+	n = t->root;
+
+	if (n->count == 1 && n->level > 0) {
+	    for (i = 0; i < t->nodecard; i++) {
+		if (n->branch[i].child)
 		    break;
 	    }
-	    assert(tmp_nptr);
-	    RTreeFreeNode(*nn);
-	    *nn = tmp_nptr;
+	    t->root = n->branch[i].child;
+	    RTreeFreeNode(n);
 	}
 	return 0;
     }
@@ -322,3 +503,20 @@ int RTreeDeleteRect1(struct Rect *R, struct Node *Child, struct Node **Nn)
 	return 1;
     }
 }
+
+/* 
+ * Delete a data rectangle from an index structure.
+ * Pass in a pointer to a Rect, the tid of the record, ptr RTree.
+ * Returns 1 if record not found, 0 if success.
+ * RTreeDeleteRect1 provides for eliminating the root.
+ *
+ * RTreeDeleteRect() should be called by external functions instead of 
+ * RTreeDeleteRect1()
+ * wrapper for RTreeDeleteRect1 not really needed, but restricts 
+ * compile warnings to rtree lib
+ * this way it's easier to fix if necessary? 
+ */
+int RTreeDeleteRect(struct Rect *r, int tid, struct RTree *t)
+{
+    return RTreeDeleteRect1(r, (struct Node *)tid, t);
+}

lib/vector/rtree/index.h

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -17,9 +18,11 @@
 #ifndef _INDEX_
 #define _INDEX_
 
+#include <sys/types.h>
+
 /* PGSIZE is normally the natural page size of the machine */
 #define PGSIZE	512
-#define NUMDIMS	3		/* number of dimensions */
+#define NUMDIMS	3		/* maximum number of dimensions */
 
 /* typedef float RectReal; */
 typedef double RectReal;
@@ -37,35 +40,77 @@ typedef double RectReal;
 
 #define NUMSIDES 2*NUMDIMS
 
+/* max branching factor of a node */
+/* was (PGSIZE-(2 * sizeof(int))) / sizeof(struct Branch)
+ * this is LFS dependent, not good
+ * on 32 bit without LFS this is 9.69
+ * on 32 bit with LFS and on 64 bit this is 9 */
+#define MAXCARD 10
+
+/* R*-tree: number of branches to be force-reinserted when adding a branch */
+#define FORCECARD 3
+
+/* maximum no of levels = tree depth */
+#define MAXLEVEL 20        /* 4^MAXLEVEL items are guaranteed to fit into the tree */
+
 struct Rect
 {
     RectReal boundary[NUMSIDES];	/* xmin,ymin,...,xmax,ymax,... */
 };
 
-struct Node;
+struct Node;               /* node for memory based index */
 
-struct Branch
+struct Branch              /* branch for memory based index */
 {
     struct Rect rect;
-    struct Node *child;
+    struct Node *child;    /* pointer to child node */
 };
 
-/* max branching factor of a node */
-#define MAXCARD (int)((PGSIZE-(2*sizeof(int))) / sizeof(struct Branch))
-
-struct Node
+struct Node             /* node for memory based index */
 {
-    int count;
+    int count;          /* number of branches */
     int level;			/* 0 is leaf, others positive */
     struct Branch branch[MAXCARD];
 };
 
+struct RTree
+{
+    /* RTree setup info */
+    unsigned char ndims;    /* number of dimensions */
+    unsigned char nsides;   /* number of sides = 2 * ndims */
+    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 */
+    unsigned int n_nodes;   /* number of nodes */
+    unsigned int n_leafs;   /* number of data items (level 0 leafs) */
+    int n_levels;           /* n levels = root level */
+    
+    /* settings for RTree building */
+    int nodecard;           /* max number of childs in node */
+    int leafcard;           /* max number of childs in leaf */
+    int min_node_fill;      /* balance criteria for node splitting */
+    int min_leaf_fill;      /* balance criteria for leaf splitting */
+    
+    struct Node *root;    /* pointer to root node */
+
+    off_t rootpos;         /* root node position in file */
+};
+
 struct ListNode
 {
     struct ListNode *next;
     struct Node *node;
 };
 
+struct ListBranch
+{
+    struct ListBranch *next;
+    struct Branch b;
+    int level;
+};
+
 /*
  * If passed to a tree search, this callback function will be called
  * with the ID of each data rect that overlaps the search rect
@@ -75,37 +120,41 @@ struct ListNode
  */
 typedef int (*SearchHitCallback) (int id, void *arg);
 
-
-extern int RTreeSearch(struct Node *, struct Rect *, SearchHitCallback,
+/* index.c */
+extern int RTreeSearch(struct RTree *, struct Rect *, SearchHitCallback,
 		       void *);
-extern int RTreeInsertRect(struct Rect *, int, struct Node **, int depth);
-extern int RTreeInsertRect1(struct Rect *, struct Node *, struct Node **, int depth);
-extern int RTreeDeleteRect(struct Rect *, int, struct Node **);
-extern int RTreeDeleteRect1(struct Rect *, struct Node *, struct Node **);
-extern struct Node *RTreeNewIndex(void);
-extern struct Node *RTreeNewNode(void);
+extern int RTreeInsertRect(struct Rect *, int, struct RTree *t);
+extern int RTreeDeleteRect(struct Rect *, int, struct RTree *t);
+extern struct RTree *RTreeNewIndex(int);
+void RTreeFreeIndex(struct RTree *);
+/* node.c */
+extern struct Node *RTreeNewNode(struct RTree *);
 extern void RTreeInitNode(struct Node *);
 extern void RTreeFreeNode(struct Node *);
-extern void RTreeDestroyNode(struct Node *);
-extern void RTreePrintNode(struct Node *, int);
+extern void RTreeDestroyNode(struct Node *, int);
+extern struct Rect RTreeNodeCover(struct Node *, struct RTree *);
+extern int RTreeAddBranch(struct Branch *, struct Node *, struct Node **, 
+            struct ListBranch **, struct Rect *, int *, struct RTree *);
+extern int RTreePickBranch(struct Rect *, struct Node *, struct RTree *);
+extern void RTreeDisconnectBranch(struct Node *, int, struct RTree *);
+extern void RTreePrintNode(struct Node *, int, struct RTree *);
 extern void RTreeTabIn(int);
-extern struct Rect RTreeNodeCover(struct Node *);
+/* rect.c */
 extern void RTreeInitRect(struct Rect *);
 extern struct Rect RTreeNullRect(void);
-extern RectReal RTreeRectArea(struct Rect *);
-extern RectReal RTreeRectSphericalVolume(struct Rect *R);
-extern RectReal RTreeRectVolume(struct Rect *R);
-extern struct Rect RTreeCombineRect(struct Rect *, struct Rect *);
-extern int RTreeOverlap(struct Rect *, struct Rect *);
+extern RectReal RTreeRectArea(struct Rect *, struct RTree *);
+extern RectReal RTreeRectSphericalVolume(struct Rect *, struct RTree *);
+extern RectReal RTreeRectVolume(struct Rect *, struct RTree *);
+extern RectReal RTreeRectMargin(struct Rect *, struct RTree *);
+extern struct Rect RTreeCombineRect(struct Rect *, struct Rect *, struct RTree *);
+extern int RTreeOverlap(struct Rect *, struct Rect *, struct RTree *);
 extern void RTreePrintRect(struct Rect *, int);
-extern int RTreeAddBranch(struct Branch *, struct Node *, struct Node **);
-extern int RTreePickBranch(struct Rect *, struct Node *);
-extern void RTreeDisconnectBranch(struct Node *, int);
-extern void RTreeSplitNode(struct Node *, struct Branch *, struct Node **);
-
-extern int RTreeSetNodeMax(int);
-extern int RTreeSetLeafMax(int);
-extern int RTreeGetNodeMax(void);
-extern int RTreeGetLeafMax(void);
+/* split.c */
+extern void RTreeSplitNode(struct Node *, struct Branch *, struct Node *, struct RTree *);
+/* card.c */
+extern int RTreeSetNodeMax(int, struct RTree *);
+extern int RTreeSetLeafMax(int, struct RTree *);
+extern int RTreeGetNodeMax(struct RTree *);
+extern int RTreeGetLeafMax(struct RTree *);
 
 #endif /* _INDEX_ */

lib/vector/rtree/node.c

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -17,22 +18,30 @@
 
 #include <stdio.h>
 #include <stdlib.h>
-#include "assert.h"
+#include <assert.h>
 #include "index.h"
 #include "card.h"
 
+
+/* rectangle distances for forced removal */
+struct dist
+{
+    int id;			/* branch id */
+    RectReal distance;		/* distance to node center */
+};
+
 /* Initialize one branch cell in a node. */
 static void RTreeInitBranch(struct Branch *b)
 {
     RTreeInitRect(&(b->rect));
     b->child = NULL;
+
 }
 
 /* Initialize a Node structure. */
-void RTreeInitNode(struct Node *N)
+void RTreeInitNode(struct Node *n)
 {
-    register struct Node *n = N;
-    register int i;
+    int i;
 
     n->count = 0;
     n->level = -1;
@@ -41,119 +50,141 @@ void RTreeInitNode(struct Node *N)
 }
 
 /* Make a new node and initialize to have all branch cells empty. */
-struct Node *RTreeNewNode(void)
+struct Node *RTreeNewNode(struct RTree *t)
 {
-    register struct Node *n;
+    struct Node *n;
 
-    /* n = new Node; */
-    n = (struct Node *)malloc(sizeof(struct Node));
+    n = (struct Node *)malloc((size_t) t->nodesize);
     assert(n);
     RTreeInitNode(n);
     return n;
 }
 
-void RTreeFreeNode(struct Node *p)
-{
-    assert(p);
-    /* delete p; */
-    free(p);
-}
-
-static void RTreePrintBranch(struct Branch *b, int depth)
-{
-    RTreePrintRect(&(b->rect), depth);
-    RTreePrintNode(b->child, depth);
-}
-
-extern void RTreeTabIn(int depth)
-{
-    int i;
-
-    for (i = 0; i < depth; i++)
-	putchar('\t');
-}
-
-/* Print out the data in a node. */
-void RTreePrintNode(struct Node *n, int depth)
+void RTreeFreeNode(struct Node *n)
 {
-    int i;
-
     assert(n);
-
-    RTreeTabIn(depth);
-    fprintf(stdout, "node");
-    if (n->level == 0)
-	fprintf(stdout, " LEAF");
-    else if (n->level > 0)
-	fprintf(stdout, " NONLEAF");
-    else
-	fprintf(stdout, " TYPE=?");
-    fprintf(stdout, "  level=%d  count=%d  address=%o\n", n->level, n->count,
-	    (unsigned)n);
-
-    for (i = 0; i < n->count; i++) {
-	if (n->level == 0) {
-	    /* RTreeTabIn(depth); */
-	    /* fprintf (stdout, "\t%d: data = %d\n", i, n->branch[i].child); */
-	}
-	else {
-	    RTreeTabIn(depth);
-	    fprintf(stdout, "branch %d\n", i);
-	    RTreePrintBranch(&n->branch[i], depth + 1);
-	}
-    }
+    free(n);
 }
 
 /*
  * Find the smallest rectangle that includes all rectangles in
  * branches of a node.
  */
-struct Rect RTreeNodeCover(struct Node *N)
+struct Rect RTreeNodeCover(struct Node *n, struct RTree *t)
 {
-    register struct Node *n = N;
-    register int i, first_time = 1;
+    int i, first_time = 1, maxkids;
     struct Rect r;
 
     assert(n);
 
+    maxkids = ((n)->level > 0 ? t->nodecard : t->leafcard);
+
     RTreeInitRect(&r);
-    for (i = 0; i < MAXKIDS(n); i++)
+    for (i = 0; i < maxkids; i++)
 	if (n->branch[i].child) {
 	    if (first_time) {
 		r = n->branch[i].rect;
 		first_time = 0;
 	    }
 	    else
-		r = RTreeCombineRect(&r, &(n->branch[i].rect));
+		r = RTreeCombineRect(&r, &(n->branch[i].rect), t);
 	}
     return r;
 }
 
 /*
+ * Idea from R*-tree, modified: not overlap size but overlap number
+ * 
+ * Pick a branch from leaf nodes (current node has level 1).  Pick the 
+ * one that will result in the smallest number of overlapping siblings.
+ * This will result in the least ambiguous node covering the new 
+ * rectangle, improving search speed.
+ * In case of a tie, pick the one which needs the smallest increase in
+ * area to accomodate the new rectangle, then the smallest area before,
+ * to get the best resolution when searching.
+ */
+
+static int RTreePickBranch1(struct Rect *r, struct Node *n, struct RTree *t)
+{
+    struct Rect *rr;
+    int i, j;
+    RectReal increase, bestIncr = (RectReal) - 1, area, bestArea = 0;
+    int best = 0, bestoverlap;
+    struct Rect tmp_rect;
+    int overlap;
+
+    assert(r && n && t);
+
+    bestoverlap = t->nodecard + 1;
+
+    /* get the branch that will overlap with the smallest number of 
+     * sibling branches when including the new rectangle */
+    for (i = 0; i < t->nodecard; i++) {
+	if (n->branch[i].child) {
+	    rr = &n->branch[i].rect;
+	    tmp_rect = RTreeCombineRect(r, rr, t);
+	    area = RTreeRectSphericalVolume(rr, t);
+	    increase = RTreeRectSphericalVolume(&tmp_rect, 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);
+		}
+	    }
+
+	    if (overlap < bestoverlap) {
+		best = i;
+		bestoverlap = overlap;
+		bestArea = area;
+		bestIncr = increase;
+	    }
+	    else if (overlap == bestoverlap) {
+		/* resolve ties */
+		if (increase < bestIncr) {
+		    best = i;
+		    bestArea = area;
+		    bestIncr = increase;
+		}
+		else if (increase == bestIncr && area < bestArea) {
+		    best = i;
+		    bestArea = area;
+		}
+	    }
+	}
+    }
+
+    return best;
+}
+
+/*
  * Pick a branch.  Pick the one that will need the smallest increase
  * in area to accomodate the new rectangle.  This will result in the
  * least total area for the covering rectangles in the current node.
  * In case of a tie, pick the one which was smaller before, to get
  * the best resolution when searching.
  */
-int RTreePickBranch(struct Rect *R, struct Node *N)
+int RTreePickBranch(struct Rect *r, struct Node *n, struct RTree *t)
 {
-    register struct Rect *r = R;
-    register struct Node *n = N;
-    register struct Rect *rr;
-    register int i, first_time = 1;
+    struct Rect *rr;
+    int i, first_time = 1;
     RectReal increase, bestIncr = (RectReal) - 1, area, bestArea = 0;
     int best = 0;
     struct Rect tmp_rect;
 
-    assert(r && n);
+    assert(r && n && t);
+    assert((n)->level > 0);	/* must not be called on leaf node */
 
-    for (i = 0; i < MAXKIDS(n); i++) {
+    if ((n)->level == 1)
+	return RTreePickBranch1(r, n, t);
+
+    for (i = 0; i < t->nodecard; i++) {
 	if (n->branch[i].child) {
 	    rr = &n->branch[i].rect;
-	    area = RTreeRectSphericalVolume(rr);
-	    tmp_rect = RTreeCombineRect(r, rr);
-	    increase = RTreeRectSphericalVolume(&tmp_rect) - area;
+	    area = RTreeRectSphericalVolume(rr, t);
+	    tmp_rect = RTreeCombineRect(r, rr, t);
+	    increase = RTreeRectSphericalVolume(&tmp_rect, t) - area;
 	    if (increase < bestIncr || first_time) {
 		best = i;
 		bestArea = area;
@@ -163,32 +194,276 @@ int RTreePickBranch(struct Rect *R, struct Node *N)
 	    else if (increase == bestIncr && area < bestArea) {
 		best = i;
 		bestArea = area;
-		bestIncr = increase;
 	    }
 	}
     }
     return best;
 }
 
+/* Disconnect a dependent node. */
+void RTreeDisconnectBranch(struct Node *n, int i, struct RTree *t)
+{
+    int maxkids;
+
+    maxkids = ((n)->level > 0 ? t->nodecard : t->leafcard);
+
+    assert(n && i >= 0 && i < maxkids);
+    assert(n->branch[i].child);
+
+    RTreeInitBranch(&(n->branch[i]));
+    n->count--;
+}
+
+/* Destroy (free) node recursively. */
+/* NOTE: only needed for memory based index */
+void RTreeDestroyNode(struct Node *n, int nodes)
+{
+    int i;
+
+    if (n->level > 0) {		/* it is not leaf -> destroy childs */
+	for (i = 0; i < nodes; i++) {
+	    if (n->branch[i].child) {
+		RTreeDestroyNode(n->branch[i].child, nodes);
+	    }
+	}
+    }
+
+    /* Free this node */
+    RTreeFreeNode(n);
+}
+
+/**********************************************************************
+ *                                                                    *
+ *   R*-tree: force remove p (currently 3) branches for reinsertion   *
+ *                                                                    *
+ **********************************************************************/
+
+/*
+ * swap dist structs
+ */
+static void RTreeSwapDist(struct dist *a, struct dist *b)
+{
+    struct dist c;
+
+    c = *a;
+    *a = *b;
+    *b = c;
+}
+
+/*
+ * check if dist is sorted ascending to distance
+ */
+static int RTreeDistIsSorted(struct dist *d, int first, int last)
+{
+    int i;
+
+    for (i = first; i < last; i++) {
+	if (d[i].distance > d[i + 1].distance)
+	    return 0;
+    }
+
+    return 1;
+}
+
+/*
+ * partition dist for quicksort on distance
+ */
+static int RTreePartitionDist(struct dist *d, int first, int last)
+{
+    int pivot, mid = (first + last) / 2;
+    int larger, smaller;
+
+    if (last - first == 1) {	/* only two items in list */
+	if (d[first].distance > d[last].distance) {
+	    RTreeSwapDist(&(d[first]), &(d[last]));
+	}
+	return last;
+    }
+
+    /* Larger of two */
+    if (d[first].distance > d[mid].distance) {
+	larger = pivot = first;
+	smaller = mid;
+    }
+    else {
+	larger = pivot = mid;
+	smaller = first;
+    }
+
+    if (d[larger].distance > d[last].distance) {
+	/* larger is largest, get the larger of smaller and last */
+	if (d[smaller].distance > d[last].distance) {
+	    pivot = smaller;
+	}
+	else {
+	    pivot = last;
+	}
+    }
+
+    if (pivot != last) {
+	RTreeSwapDist(&(d[pivot]), &(d[last]));
+    }
+
+    pivot = first;
+
+    while (first < last) {
+	if (d[first].distance <= d[last].distance) {
+	    if (pivot != first) {
+		RTreeSwapDist(&(d[pivot]), &(d[first]));
+	    }
+	    pivot++;
+	}
+	++first;
+    }
+
+    if (pivot != last) {
+	RTreeSwapDist(&(d[pivot]), &(d[last]));
+    }
+
+    return pivot;
+}
+
+/*
+ * quicksort dist struct ascending by distance
+ * n is last valid index
+ */
+static void RTreeQuicksortDist(struct dist *d, int n)
+{
+    int pivot, first, last;
+    int s_first[MAXCARD + 1], s_last[MAXCARD + 1], stacksize;
+
+    s_first[0] = 0;
+    s_last[0] = n;
+
+    stacksize = 1;
+
+    /* use stack */
+    while (stacksize) {
+	stacksize--;
+	first = s_first[stacksize];
+	last = s_last[stacksize];
+	if (first < last) {
+	    if (!RTreeDistIsSorted(d, first, last)) {
+
+		pivot = RTreePartitionDist(d, first, last);
+
+		s_first[stacksize] = first;
+		s_last[stacksize] = pivot - 1;
+		stacksize++;
+
+		s_first[stacksize] = pivot + 1;
+		s_last[stacksize] = last;
+		stacksize++;
+	    }
+	}
+    }
+}
+
+/*
+ * Allocate space for a branch in the list used in InsertRect to
+ * store branches of nodes that are too full.
+ */
+static struct ListBranch *RTreeNewListBranch(void)
+{
+    return (struct ListBranch *)malloc(sizeof(struct ListBranch));
+    /* return new ListBranch; */
+}
+
+/*
+ * Remove branches from a node. Select the 3 branches whose rectangle 
+ * center is farthest away from node cover center.
+ * Old node updated.
+ */
+
+/* 
+ * Add a branch to the reinsertion list.  It will later
+ * be reinserted into the index structure.
+ */
+static void RTreeReInsertBranch(struct Branch b, int level,
+				struct ListBranch **ee)
+{
+    register struct ListBranch *l;
+
+    l = RTreeNewListBranch();
+    l->b = b;
+    l->level = level;
+    l->next = *ee;
+    *ee = l;
+}
+
+static void RTreeRemoveBranches(struct Node *n, struct Branch *b,
+				struct ListBranch **ee, struct Rect *cover,
+				struct RTree *t)
+{
+    int i, j, maxkids;
+    RectReal center_n[NUMDIMS], center_r, delta;
+    struct Branch branchbuf[MAXCARD + 1];
+    struct dist rdist[MAXCARD + 1];
+
+    maxkids = t->nodecard;
+
+    assert(n->count == maxkids);	/* must be full */
+
+    /* center coords of node cover */
+    for (j = 0; j < t->ndims; j++) {
+	center_n[j] = (cover->boundary[j + NUMDIMS] + 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];
+	rdist[i].distance = 0;
+	rdist[i].id = i;
+	for (j = 0; j < t->ndims; j++) {
+	    center_r =
+		(n->branch[i].rect.boundary[j + NUMDIMS] +
+		 n->branch[i].rect.boundary[j]) / 2;
+	    delta = center_n[j] - center_r;
+	    rdist[i].distance += delta * delta;
+	}
+
+	RTreeInitBranch(&(n->branch[i]));
+    }
+
+    /* new branch */
+    branchbuf[maxkids] = *b;
+    rdist[maxkids].distance = 0;
+    rdist[maxkids].id = i;
+
+    /* quicksort dist */
+    RTreeQuicksortDist(rdist, maxkids);
+
+    /* put largest three in branch list */
+    for (i = 0; i < FORCECARD; i++) {
+	RTreeReInsertBranch(branchbuf[rdist[maxkids - i].id], n->level, ee);
+    }
+    /* put remaining in node, closest to center first */
+    for (i = 0; i < maxkids - FORCECARD + 1; i++) {
+	n->branch[i] = branchbuf[rdist[i].id];
+    }
+    n->count = maxkids - FORCECARD + 1;
+}
+
 /*
  * Add a branch to a node.  Split the node if necessary.
  * Returns 0 if node not split.  Old node updated.
  * Returns 1 if node split, sets *new_node to address of new node.
  * Old node updated, becomes one of two.
  */
-int RTreeAddBranch(struct Branch *B, struct Node *N, struct Node **New_node)
+int RTreeAddBranch(struct Branch *b, struct Node *n,
+		   struct Node **new_node, struct ListBranch **ee,
+		   struct Rect *cover, int *overflow, struct RTree *t)
 {
-    register struct Branch *b = B;
-    register struct Node *n = N;
-    register struct Node **new_node = New_node;
-    register int i;
+    int i, maxkids;
 
     assert(b);
     assert(n);
 
-    if (n->count < MAXKIDS(n)) {	/* split won't be necessary */
-	for (i = 0; i < MAXKIDS(n); i++) {	/* find empty branch */
-	    if (n->branch[i].child == NULL) {
+    maxkids = ((n)->level > 0 ? t->nodecard : t->leafcard);
+
+    if (n->count < maxkids) {	/* split won't be necessary */
+	for (i = 0; i < maxkids; i++) {	/* find empty branch */
+	    if (n->branch[i].child == 0) {
 		n->branch[i] = *b;
 		n->count++;
 		break;
@@ -197,35 +472,68 @@ int RTreeAddBranch(struct Branch *B, struct Node *N, struct Node **New_node)
 	return 0;
     }
     else {
-	assert(new_node);
-	RTreeSplitNode(n, b, new_node);
-	return 1;
+	if (n->level < t->n_levels && overflow[n->level]) {
+	    /* R*-tree forced reinsert */
+	    RTreeRemoveBranches(n, b, ee, cover, t);
+	    overflow[n->level] = 0;
+	    return 2;
+	}
+	else {
+	    *new_node = RTreeNewNode(t);
+	    RTreeSplitNode(n, b, *new_node, t);
+	    return 1;
+	}
     }
 }
 
-/* Disconnect a dependent node. */
-void RTreeDisconnectBranch(struct Node *n, int i)
+/*
+ * for debugging only: print items to stdout
+ */
+
+void RTreeTabIn(int depth)
 {
-    assert(n && i >= 0 && i < MAXKIDS(n));
-    assert(n->branch[i].child);
+    int i;
 
-    RTreeInitBranch(&(n->branch[i]));
-    n->count--;
+    for (i = 0; i < depth; i++)
+	putchar('\t');
 }
 
-/* Destroy (free) node recursively. */
-void RTreeDestroyNode(struct Node *n)
+static void RTreePrintBranch(struct Branch *b, int depth, struct RTree *t)
 {
-    int i;
+    RTreePrintRect(&(b->rect), depth);
+    RTreePrintNode(b->child, depth, t);
+}
 
-    if (n->level > 0) {		/* it is not leaf -> destroy childs */
-	for (i = 0; i < NODECARD; i++) {
-	    if (n->branch[i].child) {
-		RTreeDestroyNode(n->branch[i].child);
-	    }
+/* Print out the data in a node. */
+void RTreePrintNode(struct Node *n, int depth, struct RTree *t)
+{
+    int i, maxkids;
+
+    RTreeTabIn(depth);
+
+
+    maxkids = (n->level > 0 ? t->nodecard : t->leafcard);
+
+    fprintf(stdout, "node");
+    if (n->level == 0)
+	fprintf(stdout, " LEAF");
+    else if (n->level > 0)
+	fprintf(stdout, " NONLEAF");
+    else
+	fprintf(stdout, " TYPE=?");
+    fprintf(stdout, "  level=%d  count=%d", n->level, n->count);
+
+    for (i = 0; i < maxkids; i++) {
+	if (n->level == 0) {
+	    RTreeTabIn(depth);
+	    RTreePrintRect(&(n->branch[i].rect), depth);
+	    fprintf(stdout, "\t%d: data id = %o\n", i,
+		    (unsigned)n->branch[i].child);
+	}
+	else {
+	    RTreeTabIn(depth);
+	    fprintf(stdout, "branch %d\n", i);
+	    RTreePrintBranch(&(n->branch[i]), depth + 1, t);
 	}
     }
-
-    /* Free this node */
-    RTreeFreeNode(n);
 }

lib/vector/rtree/rect.c

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -17,12 +18,11 @@
 
 #include <stdio.h>
 #include <stdlib.h>
-#include "assert.h"
+#include <assert.h>
 #include "index.h"
 
 #include <float.h>
 #include <math.h>
-#include <grass/gis.h>
 
 #define BIG_NUM (FLT_MAX/4.0)
 
@@ -145,7 +145,7 @@ void RTreePrintRect(struct Rect *R, int depth)
 /*-----------------------------------------------------------------------------
 | Calculate the n-dimensional volume of a rectangle
 -----------------------------------------------------------------------------*/
-RectReal RTreeRectVolume(struct Rect *R)
+RectReal RTreeRectVolume(struct Rect *R, struct RTree *t)
 {
     register struct Rect *r = R;
     register int i;
@@ -155,7 +155,7 @@ RectReal RTreeRectVolume(struct Rect *R)
     if (Undefined(r))
 	return (RectReal) 0;
 
-    for (i = 0; i < NUMDIMS; i++)
+    for (i = 0; i < t->ndims; i++)
 	volume *= r->boundary[i + NUMDIMS] - r->boundary[i];
     assert(volume >= 0.0);
     return volume;
@@ -229,7 +229,7 @@ const double UnitSphereVolumes[] = {
  * A fast approximation to the volume of the bounding sphere for the
  * given Rect. By Paul B.
  */
-RectReal RTreeRectSphericalVolume(struct Rect *R)
+RectReal RTreeRectSphericalVolume(struct Rect *R, struct RTree *t)
 {
     register struct Rect *r = R;
     register int i;
@@ -238,57 +238,56 @@ RectReal RTreeRectSphericalVolume(struct Rect *R)
     assert(r);
     if (Undefined(r))
 	return (RectReal) 0;
-    for (i = 0; i < NUMDIMS; i++) {
+    for (i = 0; i < t->ndims; i++) {
 	c_size = r->boundary[i + NUMDIMS] - r->boundary[i];
 	if (c_size > maxsize)
 	    maxsize = c_size;
     }
-    return (RectReal) (pow(maxsize / 2, NUMDIMS) * UnitSphereVolume);
+    return (RectReal) (pow(maxsize / 2, NUMDIMS) *
+		       UnitSphereVolumes[t->ndims]);
 }
 #endif
 
 /*
  * The exact volume of the bounding sphere for the given Rect.
  */
-RectReal RTreeRectSphericalVolume(struct Rect * R)
+RectReal RTreeRectSphericalVolume(struct Rect * r, struct RTree * t)
 {
-    register struct Rect *r = R;
-    register int i;
-    register double sum_of_squares = 0, radius;
+    int i;
+    double sum_of_squares = 0, radius, half_extent;
 
     assert(r);
     if (Undefined(r))
 	return (RectReal) 0;
-    for (i = 0; i < NUMDIMS; i++) {
-	double half_extent = (r->boundary[i + NUMDIMS] - r->boundary[i]) / 2;
+    for (i = 0; i < t->ndims; i++) {
+	half_extent = (r->boundary[i + NUMDIMS] - r->boundary[i]) / 2;
 
 	sum_of_squares += half_extent * half_extent;
     }
     radius = sqrt(sum_of_squares);
-    return (RectReal) (pow(radius, NUMDIMS) * UnitSphereVolume);
+    return (RectReal) (pow(radius, t->ndims) * UnitSphereVolumes[t->ndims]);
 }
 
 
 /*-----------------------------------------------------------------------------
 | Calculate the n-dimensional surface area of a rectangle
 -----------------------------------------------------------------------------*/
-RectReal RTreeRectSurfaceArea(struct Rect * R)
+RectReal RTreeRectSurfaceArea(struct Rect * r, struct RTree * t)
 {
-    register struct Rect *r = R;
-    register int i, j;
-    register RectReal sum = (RectReal) 0;
+    int i, j;
+    RectReal j_extent, sum = (RectReal) 0;
 
     assert(r);
     if (Undefined(r))
 	return (RectReal) 0;
 
-    for (i = 0; i < NUMDIMS; i++) {
+    for (i = 0; i < t->ndims; i++) {
 	RectReal face_area = (RectReal) 1;
 
-	for (j = 0; j < NUMDIMS; j++)
+	for (j = 0; j < t->ndims; j++)
 	    /* exclude i extent from product in this dimension */
 	    if (i != j) {
-		RectReal j_extent = r->boundary[j + NUMDIMS] - r->boundary[j];
+		j_extent = r->boundary[j + NUMDIMS] - r->boundary[j];
 
 		face_area *= j_extent;
 	    }
@@ -298,14 +297,31 @@ RectReal RTreeRectSurfaceArea(struct Rect * R)
 }
 
 
+/*-----------------------------------------------------------------------------
+| Calculate the n-dimensional margin of a rectangle
+| the margin is the sum of the lengths of the edges
+-----------------------------------------------------------------------------*/
+RectReal RTreeRectMargin(struct Rect * r, struct RTree * t)
+{
+    int i;
+    RectReal margin = 0.0;
+
+    assert(r);
+
+    for (i = 0; i < t->ndims; i++) {
+	margin += r->boundary[i + NUMDIMS] - r->boundary[i];
+    }
+
+    return margin;
+}
+
 
 /*-----------------------------------------------------------------------------
 | Combine two rectangles, make one that includes both.
 -----------------------------------------------------------------------------*/
-struct Rect RTreeCombineRect(struct Rect *R, struct Rect *Rr)
+struct Rect RTreeCombineRect(struct Rect *r, struct Rect *rr, struct RTree *t)
 {
-    register struct Rect *r = R, *rr = Rr;
-    register int i, j;
+    int i, j;
     struct Rect new_rect;
 
     assert(r && rr);
@@ -316,7 +332,7 @@ struct Rect RTreeCombineRect(struct Rect *R, struct Rect *Rr)
     if (Undefined(rr))
 	return *r;
 
-    for (i = 0; i < NUMDIMS; i++) {
+    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]);
@@ -328,14 +344,13 @@ struct Rect RTreeCombineRect(struct Rect *R, struct Rect *Rr)
 /*-----------------------------------------------------------------------------
 | Decide whether two rectangles overlap.
 -----------------------------------------------------------------------------*/
-int RTreeOverlap(struct Rect *R, struct Rect *S)
+int RTreeOverlap(struct Rect *r, struct Rect *s, struct RTree *t)
 {
-    register struct Rect *r = R, *s = S;
     register int i, j;
 
     assert(r && s);
 
-    for (i = 0; i < NUMDIMS; i++) {
+    for (i = 0; i < t->ndims; i++) {
 	j = i + NUMDIMS;	/* index for high sides */
 	if (r->boundary[i] > s->boundary[j] ||
 	    s->boundary[i] > r->boundary[j]) {
@@ -349,12 +364,11 @@ int RTreeOverlap(struct Rect *R, struct Rect *S)
 /*-----------------------------------------------------------------------------
 | Decide whether rectangle r is contained in rectangle s.
 -----------------------------------------------------------------------------*/
-int RTreeContained(struct Rect *R, struct Rect *S)
+int RTreeContained(struct Rect *r, struct Rect *s, struct RTree *t)
 {
-    register struct Rect *r = R, *s = S;
     register int i, j, result;
 
-    assert(r && s); /* same as in RTreeOverlap() */
+    assert(r && s);		/* same as in RTreeOverlap() */
 
     /* undefined rect is contained in any other */
     if (Undefined(r))
@@ -365,7 +379,7 @@ int RTreeContained(struct Rect *R, struct Rect *S)
 	return FALSE;
 
     result = TRUE;
-    for (i = 0; i < NUMDIMS; i++) {
+    for (i = 0; i < t->ndims; i++) {
 	j = i + NUMDIMS;	/* index for high sides */
 	result = result && r->boundary[i] >= s->boundary[i]
 	    && r->boundary[j] <= s->boundary[j];

lib/vector/rtree/rtree.h

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -17,5 +18,4 @@
 
 #include <grass/rtree/card.h>
 #include <grass/rtree/index.h>
-/* #include <rtree/split_l.h> */
-#include <grass/rtree/split_q.h>
+#include <grass/rtree/split.h>

lib/vector/rtree/split.c

@@ -0,0 +1,621 @@
+
+/***********************************************************************
+* MODULE:       R-Tree library 
+*              
+* AUTHOR(S):    Antonin Guttman - original code
+*               Daniel Green (green@superliminal.com) - major clean-up
+*                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
+*               
+* PURPOSE:      Multidimensional index
+*
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+*
+*               This program is free software under the GNU General 
+* 		Public License (>=v2). Read the file COPYING that comes 
+* 		with GRASS for details.
+* 
+***********************************************************************/
+
+#include <stdio.h>
+#include <assert.h>
+#include "index.h"
+#include "card.h"
+#include "split.h"
+
+struct Branch BranchBuf[MAXCARD + 1];
+int BranchCount;
+struct Rect CoverSplit;
+RectReal CoverSplitArea;
+
+/* variables for finding a partition */
+struct PartitionVars Partitions[METHODS];
+
+/*----------------------------------------------------------------------
+| Load branch buffer with branches from full node plus the extra branch.
+----------------------------------------------------------------------*/
+static void RTreeGetBranches(struct Node *n, struct Branch *b,
+			     struct RTree *t)
+{
+    int i, maxkids;
+
+    assert(n);
+    assert(b);
+
+    maxkids = ((n)->level > 0 ? t->nodecard : t->leafcard);
+
+    /* load the branch buffer */
+    for (i = 0; i < maxkids; i++) {
+	assert(n->branch[i].child);	/* n should have every entry full */
+	BranchBuf[i] = n->branch[i];
+    }
+    BranchBuf[maxkids] = *b;
+    BranchCount = maxkids + 1;
+
+    /* calculate rect containing all in the set */
+    CoverSplit = BranchBuf[0].rect;
+    for (i = 1; i < maxkids + 1; i++) {
+	CoverSplit = RTreeCombineRect(&CoverSplit, &BranchBuf[i].rect, t);
+    }
+    CoverSplitArea = RTreeRectSphericalVolume(&CoverSplit, t);
+
+    RTreeInitNode(n);
+}
+
+/*----------------------------------------------------------------------
+| Put a branch in one of the groups.
+----------------------------------------------------------------------*/
+static void RTreeClassify(int i, int group, struct PartitionVars *p,
+			  struct RTree *t)
+{
+    assert(p);
+    assert(!p->taken[i]);
+
+    p->partition[i] = group;
+    p->taken[i] = TRUE;
+
+    if (p->count[group] == 0)
+	p->cover[group] = BranchBuf[i].rect;
+    else
+	p->cover[group] =
+	    RTreeCombineRect(&BranchBuf[i].rect, &p->cover[group], t);
+    p->area[group] = RTreeRectSphericalVolume(&p->cover[group], t);
+    p->count[group]++;
+}
+
+/**********************************************************************
+ *                                                                    *
+ *            Toni Guttmann's quadratic splitting method              *
+ *                                                                    *
+ **********************************************************************/
+
+/*----------------------------------------------------------------------
+| Pick two rects from set to be the first elements of the two groups.
+| Pick the two that waste the most area if covered by a single
+| rectangle.
+----------------------------------------------------------------------*/
+static void RTreePickSeeds(struct PartitionVars *p, 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);
+
+    worst = -CoverSplitArea - 1;
+    for (i = 0; i < p->total - 1; i++) {
+	for (j = i + 1; j < p->total; j++) {
+	    struct Rect one_rect;
+
+	    one_rect = RTreeCombineRect(&BranchBuf[i].rect,
+					&BranchBuf[j].rect, t);
+	    waste =
+		RTreeRectSphericalVolume(&one_rect, t) - area[i] - area[j];
+	    if (waste > worst) {
+		worst = waste;
+		seed0 = i;
+		seed1 = j;
+	    }
+	}
+    }
+    RTreeClassify(seed0, 0, p, t);
+    RTreeClassify(seed1, 1, p, t);
+}
+
+/*----------------------------------------------------------------------
+| Copy branches from the buffer into two nodes according to the
+| partition.
+----------------------------------------------------------------------*/
+static void RTreeLoadNodes(struct Node *n, struct Node *q,
+			   struct PartitionVars *p, struct RTree *t)
+{
+    int i;
+
+    assert(n);
+    assert(q);
+    assert(p);
+    assert(t);
+
+    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);
+	else if (p->partition[i] == 1)
+	    RTreeAddBranch(&BranchBuf[i], q, NULL, NULL, NULL, NULL, t);
+    }
+}
+
+/*----------------------------------------------------------------------
+| Initialize a PartitionVars structure.
+----------------------------------------------------------------------*/
+void RTreeInitPVars(struct PartitionVars *p, int maxrects, int minfill)
+{
+    int i;
+
+    assert(p);
+
+    p->count[0] = p->count[1] = 0;
+    p->cover[0] = p->cover[1] = RTreeNullRect();
+    p->area[0] = p->area[1] = (RectReal) 0;
+    p->total = maxrects;
+    p->minfill = minfill;
+    for (i = 0; i < maxrects; i++) {
+	p->taken[i] = FALSE;
+	p->partition[i] = -1;
+    }
+}
+
+/*----------------------------------------------------------------------
+| Print out data for a partition from PartitionVars struct.
+| Unused, for debugging only
+----------------------------------------------------------------------*/
+static void RTreePrintPVars(struct PartitionVars *p)
+{
+    int i;
+
+    assert(p);
+
+    fprintf(stdout, "\npartition:\n");
+    for (i = 0; i < p->total; i++) {
+	fprintf(stdout, "%3d\t", i);
+    }
+    fprintf(stdout, "\n");
+    for (i = 0; i < p->total; i++) {
+	if (p->taken[i])
+	    fprintf(stdout, "  t\t");
+	else
+	    fprintf(stdout, "\t");
+    }
+    fprintf(stdout, "\n");
+    for (i = 0; i < p->total; i++) {
+	fprintf(stdout, "%3d\t", p->partition[i]);
+    }
+    fprintf(stdout, "\n");
+
+    fprintf(stdout, "count[0] = %d  area = %f\n", p->count[0], p->area[0]);
+    fprintf(stdout, "count[1] = %d  area = %f\n", p->count[1], p->area[1]);
+    if (p->area[0] + p->area[1] > 0) {
+	fprintf(stdout, "total area = %f  effectiveness = %3.2f\n",
+		p->area[0] + p->area[1],
+		(float)CoverSplitArea / (p->area[0] + p->area[1]));
+    }
+    fprintf(stdout, "cover[0]:\n");
+    RTreePrintRect(&p->cover[0], 0);
+
+    fprintf(stdout, "cover[1]:\n");
+    RTreePrintRect(&p->cover[1], 0);
+}
+
+/*----------------------------------------------------------------------
+| Method #0 for choosing a partition: this is Toni Guttmann's quadratic
+| split
+|
+| As the seeds for the two groups, pick the two rects that would waste
+| the most area if covered by a single rectangle, i.e. evidently the 
+| worst pair to have in the same group. Of the remaining, one at a time 
+| is chosen to be put in one of the two groups. The one chosen is the 
+| one with the greatest difference in area expansion depending on which
+| group - the rect most strongly attracted to one group and repelled
+| from the other. If one group gets too full (more would force other 
+| group to violate min fill requirement) then other group gets the rest.
+| These last are the ones that can go in either group most easily.
+----------------------------------------------------------------------*/
+static void RTreeMethodZero(struct PartitionVars *p, int minfill,
+			    struct RTree *t)
+{
+    int i;
+    RectReal biggestDiff;
+    int group, chosen = 0, betterGroup = 0;
+
+    assert(p);
+
+    RTreeInitPVars(p, BranchCount, minfill);
+    RTreePickSeeds(p, t);
+
+    while (p->count[0] + p->count[1] < p->total
+	   && p->count[0] < p->total - p->minfill
+	   && p->count[1] < p->total - p->minfill) {
+	biggestDiff = (RectReal) - 1.;
+	for (i = 0; i < p->total; i++) {
+	    if (!p->taken[i]) {
+		struct 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];
+		diff = growth1 - growth0;
+		if (diff >= 0)
+		    group = 0;
+		else {
+		    group = 1;
+		    diff = -diff;
+		}
+
+		if (diff > biggestDiff) {
+		    biggestDiff = diff;
+		    chosen = i;
+		    betterGroup = group;
+		}
+		else if (diff == biggestDiff &&
+			 p->count[group] < p->count[betterGroup]) {
+		    chosen = i;
+		    betterGroup = group;
+		}
+	    }
+	}
+	RTreeClassify(chosen, betterGroup, p, t);
+    }
+
+    /* if one group too full, put remaining rects in the other */
+    if (p->count[0] + p->count[1] < p->total) {
+	if (p->count[0] >= p->total - p->minfill)
+	    group = 1;
+	else
+	    group = 0;
+	for (i = 0; i < p->total; i++) {
+	    if (!p->taken[i])
+		RTreeClassify(i, group, p, t);
+	}
+    }
+
+    assert(p->count[0] + p->count[1] == p->total);
+    assert(p->count[0] >= p->minfill && p->count[1] >= p->minfill);
+}
+
+/**********************************************************************
+ *                                                                    *
+ *           Norbert Beckmann's R*-tree splitting method              *
+ *                                                                    *
+ **********************************************************************/
+
+/*----------------------------------------------------------------------
+| swap branches
+----------------------------------------------------------------------*/
+static void RTreeSwapBranches(struct Branch *a, struct Branch *b)
+{
+    struct Branch c;
+
+    c = *a;
+    *a = *b;
+    *b = c;
+}
+
+/*----------------------------------------------------------------------
+| compare branches for given rectangle side
+| return 1 if a > b
+| return 0 if a == b
+| return -1 if a < b
+----------------------------------------------------------------------*/
+static int RTreeCompareBranches(struct Branch *a, struct Branch *b, int side)
+{
+    if (a->rect.boundary[side] > b->rect.boundary[side])
+	return 1;
+    else if (a->rect.boundary[side] < b->rect.boundary[side])
+	return -1;
+
+    /* boundaries are equal */
+    return 0;
+}
+
+/*----------------------------------------------------------------------
+| check if BranchBuf is sorted along given axis (dimension)
+----------------------------------------------------------------------*/
+static int RTreeBranchBufIsSorted(int first, int last, int side)
+{
+    int i;
+
+    for (i = first; i < last; i++) {
+	if (RTreeCompareBranches(&(BranchBuf[i]), &(BranchBuf[i + 1]), side)
+	    == 1)
+	    return 0;
+    }
+
+    return 1;
+}
+
+/*----------------------------------------------------------------------
+| partition BranchBuf for quicksort along given axis (dimension)
+----------------------------------------------------------------------*/
+static int RTreePartitionBranchBuf(int first, int last, int side)
+{
+    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]));
+	}
+	return last;
+    }
+
+    /* Larger of two */
+    if (RTreeCompareBranches(&(BranchBuf[first]), &(BranchBuf[mid]), side) ==
+	1) {
+	larger = pivot = first;
+	smaller = mid;
+    }
+    else {
+	larger = pivot = mid;
+	smaller = first;
+    }
+
+    if (RTreeCompareBranches(&(BranchBuf[larger]), &(BranchBuf[last]), side)
+	== 1) {
+	/* larger is largest, get the larger of smaller and last */
+	if (RTreeCompareBranches
+	    (&(BranchBuf[smaller]), &(BranchBuf[last]), side) == 1) {
+	    pivot = smaller;
+	}
+	else {
+	    pivot = last;
+	}
+    }
+
+    if (pivot != last) {
+	RTreeSwapBranches(&(BranchBuf[pivot]), &(BranchBuf[last]));
+    }
+
+    pivot = first;
+
+    while (first < last) {
+	if (RTreeCompareBranches
+	    (&(BranchBuf[first]), &(BranchBuf[last]), side) != 1) {
+	    if (pivot != first) {
+		RTreeSwapBranches(&(BranchBuf[pivot]), &(BranchBuf[first]));
+	    }
+	    pivot++;
+	}
+	++first;
+    }
+
+    if (pivot != last) {
+	RTreeSwapBranches(&(BranchBuf[pivot]), &(BranchBuf[last]));
+    }
+
+    return pivot;
+}
+
+/*----------------------------------------------------------------------
+| quicksort BranchBuf along given side
+----------------------------------------------------------------------*/
+static void RTreeQuicksortBranchBuf(int side)
+{
+    int pivot, first, last;
+    int s_first[MAXCARD + 1], s_last[MAXCARD + 1], stacksize;
+
+    s_first[0] = 0;
+    s_last[0] = BranchCount - 1;
+
+    stacksize = 1;
+
+    /* use stack */
+    while (stacksize) {
+	stacksize--;
+	first = s_first[stacksize];
+	last = s_last[stacksize];
+	if (first < last) {
+	    if (!RTreeBranchBufIsSorted(first, last, side)) {
+
+		pivot = RTreePartitionBranchBuf(first, last, side);
+
+		s_first[stacksize] = first;
+		s_last[stacksize] = pivot - 1;
+		stacksize++;
+
+		s_first[stacksize] = pivot + 1;
+		s_last[stacksize] = last;
+		stacksize++;
+	    }
+	}
+    }
+}
+
+/*----------------------------------------------------------------------
+| Method #1 for choosing a partition: this is the R*-tree method
+|
+| Pick the axis with the smallest margin increase (keep rectangles 
+| square).
+| Along the chosen split axis, choose the distribution with the mimimum 
+| overlap-value. Resolve ties by choosing the distribution with the
+| mimimum area-value.
+| If one group gets too full (more would force other group to violate min
+| fill requirement) then other group gets the rest.
+| These last are the ones that can go in either group most easily.
+----------------------------------------------------------------------*/
+static void RTreeMethodOne(struct PartitionVars *p, int minfill,
+			   struct RTree *t)
+{
+    int i, j, k, l, s, maxkids, first_time = 1;
+    int axis = 0, best_axis = 0, side = 0, best_side[NUMDIMS];
+    int best_cut[NUMDIMS];
+    RectReal margin, smallest_margin = 0;
+    struct Rect *r1, *r2, testrect1, testrect2, upperrect, orect;
+    int minfill1 = minfill - 1;
+    RectReal overlap, vol, smallest_overlap = -1, smallest_vol = -1;
+
+    assert(p);
+
+    RTreeInitPVars(p, BranchCount, minfill);
+
+    maxkids = (minfill == t->min_leaf_fill ? t->leafcard : t->nodecard);
+
+    /* choose split axis */
+    /* For each dimension, sort rectangles first by lower boundary then
+     * by upper boundary. Get the smallest margin. */
+    for (i = 0; i < t->ndims; i++) {
+	axis = i;
+	best_cut[i] = 0;
+	best_side[i] = 0;
+
+	smallest_overlap = -1;
+	smallest_vol = -1;
+
+	/* first lower then upper bounds for each axis */
+	for (s = 0; s < 2; s++) {
+	    RTreeQuicksortBranchBuf(i + s * NUMDIMS);
+
+	    side = s;
+
+	    testrect1 = BranchBuf[0].rect;
+	    upperrect = BranchBuf[maxkids].rect;
+
+	    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);
+	    }
+	    r2 = &BranchBuf[maxkids - minfill1].rect;
+	    upperrect = RTreeCombineRect(&upperrect, r2, t);
+
+	    /* check distributions for this axis, adhere the minimum node fill */
+	    for (j = minfill1; j < BranchCount - minfill; j++) {
+
+		r1 = &BranchBuf[j].rect;
+		testrect1 = RTreeCombineRect(&testrect1, r1, t);
+
+		testrect2 = upperrect;
+		for (k = j + 1; k < BranchCount - minfill; k++) {
+		    r2 = &BranchBuf[k].rect;
+		    testrect2 = RTreeCombineRect(&testrect2, r2, t);
+		}
+
+		/* the margin is the sum of the lengths of the edges of a rectangle */
+		margin =
+		    RTreeRectMargin(&testrect1,
+				    t) + RTreeRectMargin(&testrect2, t);
+
+		/* remember best axis */
+		if (margin <= smallest_margin || first_time) {
+		    smallest_margin = margin;
+		    best_axis = i;
+		    first_time = 0;
+		}
+
+		/* remember best distribution for this axis */
+
+		/* overlap size */
+		overlap = 1;
+
+		for (k = 0; k < t->ndims; k++) {
+		    /* no overlap */
+		    if (testrect1.boundary[k] >
+			testrect2.boundary[k + NUMDIMS] ||
+			testrect1.boundary[k + NUMDIMS] <
+			testrect2.boundary[k]) {
+			overlap = 0;
+			break;
+		    }
+		    /* get overlap */
+		    else {
+			if (testrect1.boundary[k] > testrect2.boundary[k])
+			    orect.boundary[k] = testrect1.boundary[k];
+			else
+			    orect.boundary[k] = testrect2.boundary[k];
+
+			l = k + NUMDIMS;
+			if (testrect1.boundary[l] < testrect2.boundary[l])
+			    orect.boundary[l] = testrect1.boundary[l];
+			else
+			    orect.boundary[l] = testrect2.boundary[l];
+		    }
+		}
+		if (overlap)
+		    overlap = RTreeRectVolume(&orect, t);
+
+		vol =
+		    RTreeRectVolume(&testrect1,
+				    t) + RTreeRectVolume(&testrect2, t);
+
+		/* get best cut for this axis */
+		if (overlap <= smallest_overlap || smallest_overlap < 0) {
+		    smallest_overlap = overlap;
+		    smallest_vol = vol;
+		    best_cut[i] = j;
+		    best_side[i] = s;
+		}
+		else if (overlap == smallest_overlap) {
+		    /* resolve ties by minimum volume */
+		    if (vol <= smallest_vol || smallest_vol < 0) {
+			smallest_vol = vol;
+			best_cut[i] = j;
+			best_side[i] = s;
+		    }
+		}
+	    }			/* end of distribution check */
+	}			/* end of side check */
+    }				/* end of axis check */
+
+    /* Use best distribution to classify branches */
+    if (best_axis != axis || best_side[best_axis] != side)
+	RTreeQuicksortBranchBuf(best_axis + best_side[best_axis] * NUMDIMS);
+
+    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++)
+	RTreeClassify(i, 1, p, t);
+
+    assert(p->count[0] + p->count[1] == p->total);
+    assert(p->count[0] >= p->minfill && p->count[1] >= p->minfill);
+}
+
+/*----------------------------------------------------------------------
+| Split a node.
+| Divides the nodes branches and the extra one between two nodes.
+| Old node is one of the new ones, and one really new one is created.
+| May use quadratic split or R*-tree split.
+----------------------------------------------------------------------*/
+void RTreeSplitNode(struct Node *n, struct Branch *b, struct Node *nn,
+		    struct RTree *t)
+{
+    struct PartitionVars *p;
+    int level;
+
+    assert(n);
+    assert(b);
+
+    /* load all the branches into a buffer, initialize old node */
+    level = n->level;
+    RTreeGetBranches(n, b, t);
+
+    /* find partition */
+    p = &Partitions[0];
+    /* Note: can't use MINFILL(n) below since n was cleared by GetBranches() */
+    /* RTreeMethodZero(p, level > 0 ? t->min_node_fill : t->min_leaf_fill, t); */
+    RTreeMethodOne(p, level > 0 ? t->min_node_fill : t->min_leaf_fill, t);
+
+    /*
+     * put branches from buffer into 2 nodes
+     * according to chosen partition
+     */
+    (nn)->level = n->level = level;
+    RTreeLoadNodes(n, nn, p, t);
+    assert(n->count + nn->count == p->total);
+}

lib/vector/rtree/split_q.h

@@ -5,10 +5,11 @@
 * AUTHOR(S):    Antonin Guttman - original code
 *               Daniel Green (green@superliminal.com) - major clean-up
 *                               and implementation of bounding spheres
+*               Markus Metz - R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
+* COPYRIGHT:    (C) 2009 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -37,3 +38,4 @@ extern RectReal CoverSplitArea;
 
 /* variables for finding a partition */
 extern struct PartitionVars Partitions[METHODS];
+

lib/vector/rtree/split_q.c

@@ -1,317 +0,0 @@
-
-/****************************************************************************
-* MODULE:       R-Tree library 
-*              
-* AUTHOR(S):    Antonin Guttman - original code
-*               Daniel Green (green@superliminal.com) - major clean-up
-*                               and implementation of bounding spheres
-*               
-* PURPOSE:      Multidimensional index
-*
-* COPYRIGHT:    (C) 2001 by the GRASS Development Team
-*
-*               This program is free software under the GNU General Public
-*               License (>=v2). Read the file COPYING that comes with GRASS
-*               for details.
-*****************************************************************************/
-
-#include <stdio.h>
-#include "assert.h"
-#include "index.h"
-#include "card.h"
-#include "split_q.h"
-
-struct Branch BranchBuf[MAXCARD + 1];
-int BranchCount;
-struct Rect CoverSplit;
-RectReal CoverSplitArea;
-
-/* variables for finding a partition */
-struct PartitionVars Partitions[METHODS];
-
-/*-----------------------------------------------------------------------------
-| Load branch buffer with branches from full node plus the extra branch.
------------------------------------------------------------------------------*/
-static void RTreeGetBranches(struct Node *n, struct Branch *b)
-{
-    int i;
-
-    assert(n);
-    assert(b);
-
-    /* load the branch buffer */
-    for (i = 0; i < MAXKIDS(n); i++) {
-	assert(n->branch[i].child);	/* n should have every entry full */
-	BranchBuf[i] = n->branch[i];
-    }
-    BranchBuf[MAXKIDS(n)] = *b;
-    BranchCount = MAXKIDS(n) + 1;
-
-    /* calculate rect containing all in the set */
-    CoverSplit = BranchBuf[0].rect;
-    for (i = 1; i < MAXKIDS(n) + 1; i++) {
-	CoverSplit = RTreeCombineRect(&CoverSplit, &BranchBuf[i].rect);
-    }
-    CoverSplitArea = RTreeRectSphericalVolume(&CoverSplit);
-
-    RTreeInitNode(n);
-}
-
-
-
-
-/*-----------------------------------------------------------------------------
-| Put a branch in one of the groups.
------------------------------------------------------------------------------*/
-static void RTreeClassify(int i, int group, struct PartitionVars *p)
-{
-    assert(p);
-    assert(!p->taken[i]);
-
-    p->partition[i] = group;
-    p->taken[i] = TRUE;
-
-    if (p->count[group] == 0)
-	p->cover[group] = BranchBuf[i].rect;
-    else
-	p->cover[group] =
-	    RTreeCombineRect(&BranchBuf[i].rect, &p->cover[group]);
-    p->area[group] = RTreeRectSphericalVolume(&p->cover[group]);
-    p->count[group]++;
-}
-
-
-
-
-/*-----------------------------------------------------------------------------
-| Pick two rects from set to be the first elements of the two groups.
-| Pick the two that waste the most area if covered by a single rectangle.
------------------------------------------------------------------------------*/
-static void RTreePickSeeds(struct PartitionVars *p)
-{
-    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);
-
-    worst = -CoverSplitArea - 1;
-    for (i = 0; i < p->total - 1; i++) {
-	for (j = i + 1; j < p->total; j++) {
-	    struct Rect one_rect;
-
-	    one_rect = RTreeCombineRect(&BranchBuf[i].rect,
-					&BranchBuf[j].rect);
-	    waste = RTreeRectSphericalVolume(&one_rect) - area[i] - area[j];
-	    if (waste > worst) {
-		worst = waste;
-		seed0 = i;
-		seed1 = j;
-	    }
-	}
-    }
-    RTreeClassify(seed0, 0, p);
-    RTreeClassify(seed1, 1, p);
-}
-
-
-
-
-/*-----------------------------------------------------------------------------
-| Copy branches from the buffer into two nodes according to the partition.
------------------------------------------------------------------------------*/
-static void RTreeLoadNodes(struct Node *n, struct Node *q,
-			   struct PartitionVars *p)
-{
-    int i;
-
-    assert(n);
-    assert(q);
-    assert(p);
-
-    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);
-	else if (p->partition[i] == 1)
-	    RTreeAddBranch(&BranchBuf[i], q, NULL);
-    }
-}
-
-
-
-
-/*-----------------------------------------------------------------------------
-| Initialize a PartitionVars structure.
------------------------------------------------------------------------------*/
-static void RTreeInitPVars(struct PartitionVars *p, int maxrects, int minfill)
-{
-    int i;
-
-    assert(p);
-
-    p->count[0] = p->count[1] = 0;
-    p->cover[0] = p->cover[1] = RTreeNullRect();
-    p->area[0] = p->area[1] = (RectReal) 0;
-    p->total = maxrects;
-    p->minfill = minfill;
-    for (i = 0; i < maxrects; i++) {
-	p->taken[i] = FALSE;
-	p->partition[i] = -1;
-    }
-}
-
-
-
-
-/*-----------------------------------------------------------------------------
-| Print out data for a partition from PartitionVars struct.
------------------------------------------------------------------------------*/
-static void RTreePrintPVars(struct PartitionVars *p)
-{
-    int i;
-
-    assert(p);
-
-    fprintf(stdout, "\npartition:\n");
-    for (i = 0; i < p->total; i++) {
-	fprintf(stdout, "%3d\t", i);
-    }
-    fprintf(stdout, "\n");
-    for (i = 0; i < p->total; i++) {
-	if (p->taken[i])
-	    fprintf(stdout, "  t\t");
-	else
-	    fprintf(stdout, "\t");
-    }
-    fprintf(stdout, "\n");
-    for (i = 0; i < p->total; i++) {
-	fprintf(stdout, "%3d\t", p->partition[i]);
-    }
-    fprintf(stdout, "\n");
-
-    fprintf(stdout, "count[0] = %d  area = %f\n", p->count[0], p->area[0]);
-    fprintf(stdout, "count[1] = %d  area = %f\n", p->count[1], p->area[1]);
-    if (p->area[0] + p->area[1] > 0) {
-	fprintf(stdout, "total area = %f  effectiveness = %3.2f\n",
-		p->area[0] + p->area[1],
-		(float)CoverSplitArea / (p->area[0] + p->area[1]));
-    }
-    fprintf(stdout, "cover[0]:\n");
-    RTreePrintRect(&p->cover[0], 0);
-
-    fprintf(stdout, "cover[1]:\n");
-    RTreePrintRect(&p->cover[1], 0);
-}
-
-
-/*-----------------------------------------------------------------------------
-| Method #0 for choosing a partition:
-| As the seeds for the two groups, pick the two rects that would waste the
-| most area if covered by a single rectangle, i.e. evidently the worst pair
-| to have in the same group.
-| Of the remaining, one at a time is chosen to be put in one of the two groups.
-| The one chosen is the one with the greatest difference in area expansion
-| depending on which group - the rect most strongly attracted to one group
-| and repelled from the other.
-| If one group gets too full (more would force other group to violate min
-| fill requirement) then other group gets the rest.
-| These last are the ones that can go in either group most easily.
------------------------------------------------------------------------------*/
-static void RTreeMethodZero(struct PartitionVars *p, int minfill)
-{
-    int i;
-    RectReal biggestDiff;
-    int group, chosen = 0, betterGroup = 0;
-
-    assert(p);
-
-    RTreeInitPVars(p, BranchCount, minfill);
-    RTreePickSeeds(p);
-
-    while (p->count[0] + p->count[1] < p->total
-	   && p->count[0] < p->total - p->minfill
-	   && p->count[1] < p->total - p->minfill) {
-	biggestDiff = (RectReal) - 1.;
-	for (i = 0; i < p->total; i++) {
-	    if (!p->taken[i]) {
-		struct Rect *r, rect_0, rect_1;
-		RectReal growth0, growth1, diff;
-
-		r = &BranchBuf[i].rect;
-		rect_0 = RTreeCombineRect(r, &p->cover[0]);
-		rect_1 = RTreeCombineRect(r, &p->cover[1]);
-		growth0 = RTreeRectSphericalVolume(&rect_0) - p->area[0];
-		growth1 = RTreeRectSphericalVolume(&rect_1) - p->area[1];
-		diff = growth1 - growth0;
-		if (diff >= 0)
-		    group = 0;
-		else {
-		    group = 1;
-		    diff = -diff;
-		}
-
-		if (diff > biggestDiff) {
-		    biggestDiff = diff;
-		    chosen = i;
-		    betterGroup = group;
-		}
-		else if (diff == biggestDiff &&
-			 p->count[group] < p->count[betterGroup]) {
-		    chosen = i;
-		    betterGroup = group;
-		}
-	    }
-	}
-	RTreeClassify(chosen, betterGroup, p);
-    }
-
-    /* if one group too full, put remaining rects in the other */
-    if (p->count[0] + p->count[1] < p->total) {
-	if (p->count[0] >= p->total - p->minfill)
-	    group = 1;
-	else
-	    group = 0;
-	for (i = 0; i < p->total; i++) {
-	    if (!p->taken[i])
-		RTreeClassify(i, group, p);
-	}
-    }
-
-    assert(p->count[0] + p->count[1] == p->total);
-    assert(p->count[0] >= p->minfill && p->count[1] >= p->minfill);
-}
-
-
-/*-----------------------------------------------------------------------------
-| Split a node.
-| Divides the nodes branches and the extra one between two nodes.
-| Old node is one of the new ones, and one really new one is created.
-| Tries more than one method for choosing a partition, uses best result.
------------------------------------------------------------------------------*/
-void RTreeSplitNode(struct Node *n, struct Branch *b, struct Node **nn)
-{
-    struct PartitionVars *p;
-    int level;
-
-    assert(n);
-    assert(b);
-
-    /* load all the branches into a buffer, initialize old node */
-    level = n->level;
-    RTreeGetBranches(n, b);
-
-    /* find partition */
-    p = &Partitions[0];
-    /* Note: can't use MINFILL(n) below since n was cleared by GetBranches() */
-    RTreeMethodZero(p, level > 0 ? MinNodeFill : MinLeafFill);
-
-    /*
-     * put branches from buffer into 2 nodes
-     * according to chosen partition
-     */
-    *nn = RTreeNewNode();
-    (*nn)->level = n->level = level;
-    RTreeLoadNodes(n, *nn, p);
-    assert(n->count + (*nn)->count == p->total);
-}