optionally file-based spatial index

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

lib/vector/rtree/card.c

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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/card.h

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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
@@ -23,11 +23,13 @@ extern int NODECARD;
 extern int LEAFCARD;
 
 /* balance criteria for node splitting */
-/* NOTE: can be changed if needed. */
-#define MinNodeFill (NODECARD / 2)
-#define MinLeafFill (LEAFCARD / 2)
+/* NOTE: can be changed if needed but
+ * must be >= 2 and <= (t)->[nodecard|leafcard] / 2 */
+#define MinNodeFill(t) ((t)->minfill_node_split)
+#define MinLeafFill(t) ((t)->minfill_leaf_split)
 
-#define MAXKIDS(n) ((n)->level > 0 ? NODECARD : LEAFCARD)
-#define MINFILL(n) ((n)->level > 0 ? MinNodeFill : MinLeafFill)
+#define MAXKIDS(level, t) ((level) > 0 ? (t)->nodecard : (t)->leafcard)
+#define MINFILL(level, t) ((level) > 0 ? (t)->minfill_node_split : (t)->minfill_leaf_split)
+//#define MINFILL(level, t) ((level) > 0 ? (t)->nodecard * 0.4 : (t)->leafcard * 0.4)
 
 #endif

lib/vector/rtree/gammavol.c

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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/index.c

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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,46 +17,96 @@
 *****************************************************************************/
 
 #include <stdlib.h>
+#include <sys/types.h>
+#include <unistd.h>
 #include <assert.h>
+#include <grass/gis.h>
 #include "index.h"
 #include "card.h"
 
-/* stack used for non-recursive insertion/deletion */
-struct stack
-{
-    struct Node *sn;		/* stack node */
-    int branch_id;		/* branch no to follow down */
-};
-
 /* 
  * Make a new index, empty.
+ * fp pointer to file holding index, file must be opened as w+
+ * rootpos postion of rootnode (past any header info)
  * ndims number of dimensions
  * returns pointer to RTree structure
  */
-struct RTree *RTreeNewIndex(int ndims)
+struct RTree *RTreeNewIndex(int fd, off_t rootpos, int ndims)
 {
     struct RTree *new_rtree;
     struct Node *n;
-
+    int i;
+    
     new_rtree = (struct RTree *)malloc(sizeof(struct RTree));
 
+    new_rtree->fd = fd;
+    new_rtree->rootpos = rootpos;
     new_rtree->ndims = ndims;
     new_rtree->nsides = 2 * ndims;
+    
+    new_rtree->rectsize = sizeof(struct Rect);
+
+    /* 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 Node);
     new_rtree->branchsize = sizeof(struct Branch);
-    new_rtree->rectsize = sizeof(struct Rect);
+    
+    /* create empty root node */
+    n = RTreeNewNode(new_rtree, 0);
+    new_rtree->rootlevel = n->level = 0;       /* leaf */
+    new_rtree->root = NULL;
+
+    if (fd > -1) {  /* file based */
+	/* nodecard and leafcard can be adjusted, must NOT be larger than MAXCARD */
+	new_rtree->nodecard = MAXCARD;
+	new_rtree->leafcard = MAXCARD;
+
+	/* initialize node buffer */
+	for (i = 0; i < MAXLEVEL; i++) {
+	    new_rtree->nb[i][0].dirty = 0;
+	    new_rtree->nb[i][1].dirty = 0;
+	    new_rtree->nb[i][0].pos = -1;
+	    new_rtree->nb[i][1].pos = -1;
+	    new_rtree->mru[i] = 0;
+	}
 
-    /* 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;
+	/* write empty root node */
+	lseek(new_rtree->fd, rootpos, SEEK_SET);
+	RTreeWriteNode(n, new_rtree);
+	new_rtree->nb[0][0].n = *n;
+	new_rtree->nb[0][0].pos = rootpos;
+	new_rtree->mru[0] = 0;
+	RTreeFreeNode(n);
+
+	new_rtree->insert_rect = RTreeInsertRectF;
+	new_rtree->delete_rect = RTreeDeleteRectF;
+	new_rtree->search_rect = RTreeSearchF;
+	new_rtree->valid_child = RTreeValidChildF;
+	
+    }
+    else {    /* memory based */
+	new_rtree->nodecard = MAXCARD;
+	new_rtree->leafcard = MAXCARD;
 
-    n = RTreeNewNode(new_rtree, 0);
-    new_rtree->n_levels = n->level = 0;	/* leaf */
-    new_rtree->root = n;
+	new_rtree->insert_rect = RTreeInsertRectM;
+	new_rtree->delete_rect = RTreeDeleteRectM;
+	new_rtree->search_rect = RTreeSearchM;
+	new_rtree->valid_child = RTreeValidChildM;
+
+	new_rtree->root = n;
+    }
+
+    /* minimum number of remaining children for RTreeDeleteRect */
+    /* NOTE: min fill can be changed if needed, must be < nodecard and leafcard. */
+    new_rtree->min_node_fill = (new_rtree->nodecard - 2) / 2;
+    new_rtree->min_leaf_fill = (new_rtree->leafcard - 2) / 2;
+
+    /* balance criteria for node splitting */
+    new_rtree->minfill_node_split = (new_rtree->nodecard - 1) / 2;
+    new_rtree->minfill_leaf_split = (new_rtree->leafcard - 1) / 2;
 
     new_rtree->n_nodes = 1;
     new_rtree->n_leafs = 0;
@@ -68,8 +118,12 @@ void RTreeFreeIndex(struct RTree *t)
 {
     assert(t);
 
-    if (t->root)
-	RTreeDestroyNode(t->root, t->nodecard);
+    if (t->fd > -1) {
+	if (t->free_nodes.alloc)
+	    free(t->free_nodes.pos);
+    }
+    else if (t->root)
+	RTreeDestroyNode(t->root, t->root->level ? t->nodecard : t->leafcard);
 
     free(t);
 }
@@ -79,249 +133,12 @@ void RTreeFreeIndex(struct RTree *t)
  * overlap the argument rectangle.
  * Return the number of qualifying data rects.
  */
-int RTreeSearch(struct RTree *t, struct Rect *r, SearchHitCallback shcb,
+int RTreeSearch(struct RTree *t, struct Rect *r, SearchHitCallback *shcb,
 		void *cbarg)
 {
-    struct Node *n;
-    int hitCount = 0, found;
-    int i;
-    struct stack s[MAXLEVEL];
-    int top = 0;
-
-    assert(r);
-    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.ptr &&
-		    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.ptr;
-		    s[top].branch_id = 0;
-		    found = 0;
-		    break;
-		}
-	    }
-	    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.id &&
-		    RTreeOverlap(r, &(s[top].sn->branch[i].rect), t)) {
-		    hitCount++;
-		    if (shcb) {	/* call the user-provided callback */
-			if (!shcb(s[top].sn->branch[i].child.id, 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.
- * 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, union Child child, int level,
-			    struct Node **newnode, struct RTree *t,
-			    struct ListBranch **ee, int *overflow)
-{
-    int i;
-    struct Branch b;
-    struct Node *n, *n2;
-    struct Rect *cover;
-    struct stack s[MAXLEVEL];
-    int top = 0, down = 0;
-    int result;
-
-    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.ptr;
-    }
-
-    /* Have reached level for insertion. Remove p rectangles or split */
-    if (s[top].sn->level == level) {
-	b.rect = *r;
-	/* child field of leaves contains tid of data record */
-	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);
-	/* 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.ptr = n2;
-	    b.rect = RTreeNodeCover(b.child.ptr, 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);
-
-	    /* update node count */
-	    if (result == 1) {
-		t->n_nodes++;
-	    }
-	}
-    }
-
-    *newnode = n2;
-
-    return result;
-}
-
-/* 
- * Insert a data rectangle into an index structure.
- * 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 actual insertion.
- */
-static int RTreeInsertRect1(struct Rect *r, union Child child, int level,
-			    struct RTree *t)
-{
-    struct Node *newnode;
-    struct Node *newroot;
-    struct Branch b;
-    struct ListBranch *reInsertList = NULL;
-    struct ListBranch *e;
-    int result;
-    int i, overflow[MAXLEVEL];
-
-    /* R*-tree forced reinsertion: for each level only once */
-    for (i = 0; i < MAXLEVEL; i++)
-	overflow[i] = 1;
-
-    result =
-	RTreeInsertRect2(r, child, level, &newnode, t, &reInsertList,
-			 overflow);
-			 
-    if (result == 1) {		/* root split */
-	/* grow a new root, & tree taller */
-	t->n_levels++;
-	newroot = RTreeNewNode(t, t->n_levels);
-	newroot->level = t->n_levels;
-	/* branch for old root */
-	b.rect = RTreeNodeCover(t->root, 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);
-	/* 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 */
-		t->n_levels++;
-		newroot = RTreeNewNode(t, t->n_levels);
-		newroot->level = t->n_levels;
-		/* branch for old root */
-		b.rect = RTreeNodeCover(t->root, 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);
-		/* set new root node */
-		t->root = newroot;
-		t->n_nodes++;
-	    }
-	}
-    }
+    assert(r && t);
 
-    return result;
+    return t->search_rect(t, r, shcb, cbarg);
 }
 
 /* 
@@ -333,26 +150,48 @@ static int RTreeInsertRect1(struct Rect *r, union Child child, int level,
 int RTreeInsertRect(struct Rect *r, int tid, struct RTree *t)
 {
     union Child newchild;
-    assert(r && t);
 
+    assert(r && t && tid > 0);
+    
     t->n_leafs++;
-
     newchild.id = tid;
+    
+    return t->insert_rect(r, newchild, 0, t);
+}
 
-    return RTreeInsertRect1(r, newchild, 0, t);
+/* 
+ * 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)
+{
+    union Child child;
+    
+    assert(r && t && tid > 0);
+
+    child.id = tid;
+
+    return t->delete_rect(r, child, t);
 }
 
 /*
- * Allocate space for a node in the list used in DeletRect to
+ * Allocate space for a node in the list used in DeleteRect to
  * store Nodes that are too empty.
  */
-static struct ListNode *RTreeNewListNode(void)
+struct ListNode *RTreeNewListNode(void)
 {
     return (struct ListNode *)malloc(sizeof(struct ListNode));
-    /* return new ListNode; */
 }
 
-static void RTreeFreeListNode(struct ListNode *p)
+void RTreeFreeListNode(struct ListNode *p)
 {
     free(p);
 }
@@ -361,187 +200,22 @@ static void RTreeFreeListNode(struct ListNode *p)
  * Add a node to the reinsertion list.  All its branches will later
  * be reinserted into the index structure.
  */
-static void RTreeReInsertNode(struct Node *n, struct ListNode **ee)
+void RTreeReInsertNode(struct Node *n, struct ListNode **ee)
 {
-    register struct ListNode *l;
+    struct ListNode *l = RTreeNewListNode();
 
-    l = RTreeNewListNode();
     l->node = n;
     l->next = *ee;
     *ee = l;
 }
 
-/*
- * Delete a rectangle from non-root part of an index structure.
- * Called by RTreeDeleteRect.  Descends tree non-recursively,
- * merges branches on the way back up.
- * Returns 1 if record not found, 0 if success.
+/* 
+ * Free ListBranch
  */
-static int
-RTreeDeleteRect2(struct Rect *r, union Child child, struct RTree *t,
-		 struct ListNode **ee)
+void RTreeFreeListBranch(struct ListBranch *p)
 {
-    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.ptr &&
-		    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.ptr;
-		    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.id && s[top].sn->branch[i].child.id == child.id) {	/* found item */
-		    RTreeDisconnectBranch(s[top].sn, i, t);
-		    t->n_leafs--;
-		    notfound = 0;
-		    break;
-		}
-	    }
-	    if (notfound)	/* continue searching */
-		top--;
-	}
-    }
-
-    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.ptr, ee);
-	    RTreeDisconnectBranch(s[top].sn, i, t);
-	}
-    }
-
-    return notfound;
+    free(p);
 }
 
-/*
- * 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 RTree.
- * Returns 1 if record not found, 0 if success.
- * RTreeDeleteRect1 provides for eliminating the root.
- */
-static int RTreeDeleteRect1(struct Rect *r, union Child child,
-			    struct RTree *t)
-{
-    int i, maxkids;
-    struct Node *n;
-    struct ListNode *reInsertList = NULL;
-    struct ListNode *e;
-
-    assert(r);
-    assert(t);
 
-    if (!RTreeDeleteRect2(r, child, t, &reInsertList)) {
-	/* found and deleted a data item */
-
-	/* reinsert any branches from eliminated nodes */
-	while (reInsertList) {
-	    t->n_nodes--;
-	    n = reInsertList->node;
-	    maxkids = (n->level > 0 ? t->nodecard : t->leafcard);
-	    for (i = 0; i < maxkids; i++) {
-		if (n->level > 0) {  /* reinsert node branches */
-		    if (n->branch[i].child.ptr) {
-			RTreeInsertRect1(&(n->branch[i].rect),
-					 n->branch[i].child, n->level, t);
-		    }
-		}
-		else {  /* reinsert leaf branches */
-		    if (n->branch[i].child.id) {
-			RTreeInsertRect1(&(n->branch[i].rect),
-					 n->branch[i].child, n->level, t);
-		    }
-		}
-	    }
-	    e = reInsertList;
-	    reInsertList = reInsertList->next;
-	    RTreeFreeNode(e->node);
-	    RTreeFreeListNode(e);
-	}
-
-	/* check for redundant root (not leaf, 1 child) and eliminate */
-	n = t->root;
-
-	if (n->count == 1 && n->level > 0) {
-	    for (i = 0; i < t->nodecard; i++) {
-		if (n->branch[i].child.ptr)
-		    break;
-	    }
-	    t->root = n->branch[i].child.ptr;
-	    RTreeFreeNode(n);
-	    t->n_levels--;
-	}
-	return 0;
-    }
-    else {
-	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)
-{
-    union Child child;
 
-    child.id = tid;
-    
-    return RTreeDeleteRect1(r, child, t);
-}

lib/vector/rtree/index.h

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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
@@ -18,13 +18,13 @@
 #ifndef _INDEX_
 #define _INDEX_
 
+#include <stdio.h>
 #include <sys/types.h>
 
 /* PGSIZE is normally the natural page size of the machine */
 #define PGSIZE	512
 #define NUMDIMS	3		/* maximum number of dimensions */
 
-/* typedef float RectReal; */
 typedef double RectReal;
 
 /*-----------------------------------------------------------------------------
@@ -45,43 +45,64 @@ typedef double RectReal;
  * 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 9
+#define MAXCARD 18
 
 /* R*-tree: number of branches to be force-reinserted when adding a branch */
 #define FORCECARD 2
 
 /* maximum no of levels = tree depth */
-#define MAXLEVEL 20        /* 4^MAXLEVEL items are guaranteed to fit into the tree */
+#define MAXLEVEL 20        /* 8^MAXLEVEL items are guaranteed to fit into the tree */
+
+#define NODETYPE(l, fd) ((l) == 0 ? 0 : ((fd) < 0 ? 1 : 2))
 
 struct Rect
 {
     RectReal boundary[NUMSIDES];	/* xmin,ymin,...,xmax,ymax,... */
 };
 
-struct Node;               /* node for memory based index */
+struct Node;               /* node for spatial index */
 
 union Child
 {
     int id;              /* child id */
     struct Node *ptr;    /* pointer to child node */
+    off_t pos;           /* position of child node in file */
 };
 
-struct Branch              /* branch for memory based index */
+struct Branch              /* branch for spatial index */
 {
     struct Rect rect;
     union Child child;
 };
 
-struct Node             /* node for memory based index */
+struct Node             /* node for spatial index */
 {
     int count;          /* number of branches */
-    int level;			/* 0 is leaf, others positive */
+    int level;		/* 0 is leaf, others positive */
     struct Branch branch[MAXCARD];
 };
 
+/*
+ * If passed to a tree search, this callback function will be called
+ * with the ID of each data rect that overlaps the search rect
+ * plus whatever user specific pointer was passed to the search.
+ * It can terminate the search early by returning 0 in which case
+ * the search will return the number of hits found up to that point.
+ */
+typedef int SearchHitCallback(int id, void *arg);
+
+struct RTree;
+
+typedef int rt_search_fn(struct RTree *, struct Rect *,
+                         SearchHitCallback *, void *);
+typedef int rt_insert_fn(struct Rect *, union Child, int, struct RTree *);
+typedef int rt_delete_fn(struct Rect *, union Child, struct RTree *);
+typedef int rt_valid_child_fn(union Child *);
+
 struct RTree
 {
     /* RTree setup info */
+    int fd;                 /* file descriptor */
     unsigned char ndims;    /* number of dimensions */
     unsigned char nsides;   /* number of sides = 2 * ndims */
     int nodesize;           /* node size in bytes */
@@ -89,17 +110,43 @@ struct RTree
     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 */
+    int n_nodes;            /* number of nodes */
+    int n_leafs;            /* number of data items (level 0 leafs) */
+    int rootlevel;          /* root level = tree depth */
     
     /* 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 */
+    int min_node_fill;      /* balance criteria for node removal */
+    int min_leaf_fill;      /* balance criteria for leaf removal */
+    int minfill_node_split; /* balance criteria for splitting */
+    int minfill_leaf_split; /* balance criteria for splitting */
+    
+    /* free node positions for recycling */
+    struct _recycle {
+        int avail;          /* number of available positions */
+        int alloc;          /* number of allcoated positions in *pos */
+        off_t *pos;         /* array of available positions */
+    } free_nodes;
+
+    /* node buffer for file-based index, two nodes per level
+     * more than two nodes per level would require too complex cache management:
+     * lru or pseudo-lru replacement, searching for buffered nodes */
+    struct NodeBuffer
+    {
+	struct Node n;	    /* buffered node */
+	off_t pos;	    /* file position of buffered node */
+	char dirty;         /* node in buffer was modified */
+    } nb[MAXLEVEL][2];
+    char mru[MAXLEVEL];     /* most recently used buffered node per level */
+
+    /* insert, delete, search */
+    rt_insert_fn *insert_rect;
+    rt_delete_fn *delete_rect;
+    rt_search_fn *search_rect;
+    rt_valid_child_fn *valid_child;
     
-    struct Node *root;    /* pointer to root node */
+    struct Node *root;     /* pointer to root node */
 
     off_t rootpos;         /* root node position in file */
 };
@@ -110,6 +157,12 @@ struct ListNode
     struct Node *node;
 };
 
+struct ListFNode
+{
+    struct ListFNode *next;
+    off_t node_pos;
+};
+
 struct ListBranch
 {
     struct ListBranch *next;
@@ -117,22 +170,29 @@ struct ListBranch
     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
- * plus whatever user specific pointer was passed to the search.
- * It can terminate the search early by returning 0 in which case
- * the search will return the number of hits found up to that point.
- */
-typedef int (*SearchHitCallback) (int id, void *arg);
-
 /* index.c */
-extern int RTreeSearch(struct RTree *, struct Rect *, SearchHitCallback,
+extern int RTreeSearch(struct RTree *, struct Rect *, SearchHitCallback *,
+		       void *);
+extern int RTreeInsertRect(struct Rect *, int, struct RTree *);
+extern int RTreeDeleteRect(struct Rect *, int, struct RTree *);
+extern struct RTree *RTreeNewIndex(int, off_t, int);
+extern void RTreeFreeIndex(struct RTree *);
+extern struct ListNode *RTreeNewListNode(void);
+extern void RTreeFreeListNode(struct ListNode *);
+extern void RTreeReInsertNode(struct Node *, struct ListNode **);
+extern void RTreeFreeListBranch(struct ListBranch *);
+/* indexm.c */
+extern int RTreeSearchM(struct RTree *, struct Rect *, SearchHitCallback *,
+		       void *);
+extern int RTreeInsertRectM(struct Rect *, union Child, int, struct RTree *);
+extern int RTreeDeleteRectM(struct Rect *, union Child, struct RTree *);
+extern int RTreeValidChildM(union Child *child);
+/* indexf.c */
+extern int RTreeSearchF(struct RTree *, struct Rect *, SearchHitCallback *,
 		       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 *);
+extern int RTreeInsertRectF(struct Rect *, union Child, int, struct RTree *);
+extern int RTreeDeleteRectF(struct Rect *, union Child, struct RTree *);
+extern int RTreeValidChildF(union Child *child);
 /* node.c */
 extern struct Node *RTreeNewNode(struct RTree *, int);
 extern void RTreeInitNode(struct Node *, int);
@@ -153,6 +213,7 @@ 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 RTreeCompareRect(struct Rect *, struct Rect *, struct RTree *);
 extern int RTreeOverlap(struct Rect *, struct Rect *, struct RTree *);
 extern void RTreePrintRect(struct Rect *, int);
 /* split.c */
@@ -162,5 +223,15 @@ extern int RTreeSetNodeMax(int, struct RTree *);
 extern int RTreeSetLeafMax(int, struct RTree *);
 extern int RTreeGetNodeMax(struct RTree *);
 extern int RTreeGetLeafMax(struct RTree *);
+/* fileio.c */
+extern void RTreeGetNode(struct Node *, off_t, int, struct RTree *);
+extern size_t RTreeReadNode(struct Node *, off_t, struct RTree *);
+extern void RTreePutNode(struct Node *, off_t, struct RTree *);
+extern size_t RTreeWriteNode(struct Node *, struct RTree *);
+extern size_t RTreeRewriteNode(struct Node *, off_t, struct RTree *);
+extern void RTreeUpdateRect(struct Rect *, struct Node *, off_t, int, struct RTree *);
+extern void RTreeAddNodePos(off_t, int, struct RTree *);
+extern off_t RTreeGetNodePos(struct RTree *);
+extern void RTreeFlushBuffer(struct RTree *);
 
 #endif /* _INDEX_ */

lib/vector/rtree/indexf.c

@@ -0,0 +1,472 @@
+
+/****************************************************************************
+* 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 - file-based and memory-based R*-tree
+*               
+* 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 <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <sys/types.h>
+#include <assert.h>
+#include <grass/config.h>
+#include <grass/gis.h>
+#include "index.h"
+#include "card.h"
+
+/* stack used for non-recursive insertion/deletion */
+struct fstack
+{
+    struct Node sn;	/* stack node */
+    int branch_id;	/* branch no to follow down */
+    off_t pos;		/* file position of stack node */
+};
+
+int RTreeValidChildF(union Child *child)
+{
+    return (child->pos > -1);
+}
+
+/*
+ * Search in an index tree for all data retangles that
+ * overlap the argument rectangle.
+ * Return the number of qualifying data rects.
+ */
+int RTreeSearchF(struct RTree *t, struct Rect *r,
+		 SearchHitCallback *shcb, void *cbarg)
+{
+    struct Node *n;
+    int hitCount = 0, found, currlevel;
+    int i;
+    struct fstack s[MAXLEVEL];
+    int top = 0;
+
+    /* stack size of t->rootlevel + 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 */
+    currlevel = t->rootlevel;
+    s[top].pos = t->rootpos;
+    RTreeGetNode(&(s[top].sn), s[top].pos, currlevel, t);
+    s[top].branch_id = i = 0;
+    
+    while (top >= 0) {
+	if (s[top].sn.level > 0) {		/* this is an internal node in the tree */
+	    n = &(s[top].sn);
+	    found = 1;
+	    currlevel = s[top].sn.level - 1;
+	    for (i = s[top].branch_id; i < t->nodecard; i++) {
+		if (s[top].sn.branch[i].child.pos > -1 &&
+		    RTreeOverlap(r, &(s[top].sn.branch[i].rect), t)) {
+		    s[top++].branch_id = i + 1;
+		    /* add next node to stack */
+		    s[top].pos = n->branch[i].child.pos;
+		    RTreeGetNode(&(s[top].sn), s[top].pos, currlevel, t);
+		    s[top].branch_id = 0;
+		    found = 0;
+		    break;
+		}
+	    }
+	    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.id &&
+		    RTreeOverlap(r, &(s[top].sn.branch[i].rect), t)) {
+		    hitCount++;
+		    if (shcb) {	/* call the user-provided callback */
+			if (!shcb(s[top].sn.branch[i].child.id, cbarg)) {
+			    /* callback wants to terminate search early */
+			    return hitCount;
+			}
+		    }
+		}
+	    }
+	    top--;
+	}
+    }
+
+    return hitCount;
+}
+
+/*
+ * Inserts a new data rectangle into the index structure.
+ * Non-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.
+ * 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 RTreeInsertRect2F(struct Rect *r, union Child child, int level, 
+			     struct Node *newnode, off_t *newnode_pos,
+			     struct RTree *t,
+			     struct ListBranch **ee, int *overflow)
+{
+    int i, currlevel;
+    struct Branch b;
+    struct Rect nr, *cover;
+    struct Node *n, *n2, nn;
+    struct fstack s[MAXLEVEL];
+    int top = 0, down = 0;
+    int result;
+
+    n2 = &nn;
+
+    /* add root node position to stack */
+    currlevel = t->rootlevel;
+    s[top].pos = t->rootpos;
+    RTreeGetNode(&(s[top].sn), s[top].pos, currlevel, t);
+
+    /* 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].pos = n->branch[i].child.pos;
+	currlevel--;
+	RTreeGetNode(&(s[top].sn), s[top].pos, currlevel, t);
+	assert(s[top].sn.level == currlevel);
+    }
+    
+    /* Have reached level for insertion. Add rect, split if necessary */
+    if (s[top].sn.level == level) {
+	b.rect = *r;
+	/* child field of leaves contains tid of data record */
+	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);
+	/* write out new node if node was split */
+	if (result == 1) {
+	    *newnode_pos = RTreeGetNodePos(t);
+	    RTreeWriteNode(n2, t);
+	    t->n_nodes++;
+	}
+	/* update node */
+	RTreePutNode(&(s[top].sn), s[top].pos, t);
+    }
+    else {
+	/* Not supposed to happen */
+	assert(s[top].sn.level == level);
+	return 0;
+    }
+    
+    /* go back up */
+    while (top) {
+	down = top--;
+	i = s[top].branch_id;
+	if (result == 0) {        /* branch was added */
+	    nr = RTreeCombineRect(r, &(s[top].sn.branch[i].rect), t);
+	    /* rewrite rect */
+	    if (!RTreeCompareRect(&nr, &(s[top].sn.branch[i].rect), t)) {
+		RTreeUpdateRect(&nr, &(s[top].sn), s[top].pos, i, t);
+	    }
+	}
+	else if (result == 2) {	/* branches were removed */
+	    /* get node cover of previous node */
+	    nr = RTreeNodeCover(&(s[down].sn), t);
+	    /* rewrite rect */
+	    if (!RTreeCompareRect(&nr, &(s[top].sn.branch[i].rect), t)) {
+		RTreeUpdateRect(&nr, &(s[top].sn), s[top].pos, i, 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.pos = *newnode_pos;
+	    b.rect = RTreeNodeCover(n2, 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);
+	    
+	    /* write out new node if node was split */
+	    if (result == 1) {
+		*newnode_pos = RTreeGetNodePos(t);
+		RTreeWriteNode(n2, t);
+	 	t->n_nodes++;
+	    }
+	    /* update node */
+	    RTreePutNode(&(s[top].sn), s[top].pos, t);
+	}
+    }
+
+    *newnode = *n2;
+
+    return result;
+}
+
+/* 
+ * Insert a data rectangle into an index structure.
+ * 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 actual insertion.
+ */
+int RTreeInsertRectF(struct Rect *r, union Child child, int level,
+                     struct RTree *t)
+{
+    struct Node oldroot, newroot, newnode;
+    struct Branch b;
+    struct ListBranch *e, *reInsertList = NULL;
+    int result;
+    int i, overflow[MAXLEVEL];
+    off_t newnode_pos = -1;
+
+    /* R*-tree forced reinsertion: for each level only once */
+    for (i = 0; i < MAXLEVEL; i++)
+	overflow[i] = 1;
+
+    /* no R*-tree like forced reinsertion for file-based index ? */
+    result =
+	RTreeInsertRect2F(r, child, level, &newnode, &newnode_pos, t,
+	                  &reInsertList, overflow);
+    
+    if (result == 1) {	/* root split */
+	RTreeGetNode(&oldroot, t->rootpos, t->rootlevel, t);
+	/* grow a new root, & tree taller */
+	t->rootlevel++;
+	RTreeInitNode(&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);
+	/* 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);
+	/* write new root node */
+	t->rootpos = RTreeGetNodePos(t);
+	RTreeWriteNode(&newroot, t);
+	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 =
+		RTreeInsertRect2F(&(b.rect), b.child, level, &newnode, &newnode_pos, t,
+				 &reInsertList, overflow);
+
+	    if (result == 1) {	/* root split */
+		RTreeGetNode(&oldroot, t->rootpos, t->rootlevel, t);
+		/* grow a new root, & tree taller */
+		t->rootlevel++;
+		RTreeInitNode(&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);
+		/* branch for new node created by RTreeFInsertRect2() */
+		b.rect = RTreeNodeCover(&newnode, 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);
+		t->n_nodes++;
+	    }
+	}
+    }
+
+    return result;
+}
+
+/*
+ * Delete a rectangle from non-root part of an index structure.
+ * 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
+RTreeDeleteRect2F(struct Rect *r, union Child child, struct RTree *t,
+		 struct ListNode **ee)
+{
+    int i, notfound = 1, currlevel;
+    struct Node *n;
+    struct Rect nr;
+    struct fstack s[MAXLEVEL];
+    int top = 0, down = 0;
+    int minfill;
+
+    assert(ee);
+    
+    /* add root node position to stack */
+    currlevel = t->rootlevel;
+    s[top].pos = t->rootpos;
+    RTreeGetNode(&(s[top].sn), s[top].pos, currlevel, t);
+    s[top].branch_id = 0;
+
+    while (notfound) {
+	/* go down to level 0, remember path */
+	if (s[top].sn.level > 0) {
+	    n = &(s[top].sn);
+	    currlevel = s[top].sn.level - 1;
+	    for (i = s[top].branch_id; i < t->nodecard; i++) {
+		if (n->branch[i].child.pos > -1 && RTreeOverlap(r, &(n->branch[i].rect), t)) {
+		    s[top++].branch_id = i + 1;
+		    /* add next node to stack */
+		    s[top].pos = n->branch[i].child.pos;
+		    RTreeGetNode(&(s[top].sn), s[top].pos, currlevel, t);
+		    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.id &&
+		    s[top].sn.branch[i].child.id == child.id) { /* found item */
+		    RTreeDisconnectBranch(&(s[top].sn), i, t);
+		    RTreePutNode(&(s[top].sn), s[top].pos, t);
+		    t->n_leafs--;
+		    notfound = 0;
+		    break;
+		}
+	    }
+	    if (notfound)    /* continue searching */
+		top--;
+	}
+    }
+    
+    if (notfound) {
+	return notfound;
+    }
+    
+    /* go back up */
+    while (top) {
+	down = top--;
+	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 */
+	    nr = RTreeNodeCover(&(s[down].sn), t);
+	    /* rewrite rect */
+	    if (!RTreeCompareRect(&nr, &(s[top].sn.branch[i].rect), t)) {
+		//s[top].sn.branch[i].rect = nr;
+		//RTreeRewriteRect(&nr, s[top].pos, i, t);
+		RTreeUpdateRect(&nr, &(s[top].sn), s[top].pos, i, t);
+	    }
+	}
+	else {
+	    /* 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);
+	    RTreeAddNodePos(s[top].sn.branch[i].child.pos, s[down].sn.level, t);
+	    RTreeReInsertNode(n, ee);
+	    RTreeDisconnectBranch(&(s[top].sn), i, t);
+
+	    RTreePutNode(&(s[top].sn), s[top].pos, t);
+	}
+    }
+
+    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 RTree.
+ * Returns 1 if record not found, 0 if success.
+ * RTreeDeleteRect1 provides for eliminating the root.
+ */
+int RTreeDeleteRectF(struct Rect *r, union Child child, struct RTree *t)
+{
+    int i;
+    struct Node *n, rn;
+    struct ListNode *e, *reInsertList = NULL;
+
+    if (!RTreeDeleteRect2F(r, child, t, &reInsertList)) {
+	/* found and deleted a data item */
+
+	/* reinsert any branches from eliminated nodes */
+	while (reInsertList) {
+	    t->n_nodes--;
+	    n = reInsertList->node;
+	    if (n->level > 0) {  /* reinsert node branches */
+		for (i = 0; i < t->nodecard; i++) {
+		    if (n->branch[i].child.pos > -1) {
+			RTreeInsertRectF(&(n->branch[i].rect),
+					 n->branch[i].child, n->level, t);
+		    }
+		}
+	    }
+	    else {  /* reinsert leaf branches */
+		for (i = 0; i < t->leafcard; i++) {
+		    if (n->branch[i].child.id) {
+			RTreeInsertRectF(&(n->branch[i].rect),
+					 n->branch[i].child, n->level, t);
+		    }
+		}
+	    }
+	    e = reInsertList;
+	    reInsertList = reInsertList->next;
+	    RTreeFreeNode(e->node);
+	    RTreeFreeListNode(e);
+	}
+
+	/* check for redundant root (not leaf, 1 child) and eliminate */
+	RTreeGetNode(&rn, t->rootpos, t->rootlevel, t);
+	
+	if (rn.count == 1 && rn.level > 0) {
+	    for (i = 0; i < t->nodecard; i++) {
+		if (rn.branch[i].child.pos > -1)
+		    break;
+	    }
+	    RTreeAddNodePos(t->rootpos, t->rootlevel, t);
+	    t->rootpos = rn.branch[i].child.pos;
+	    t->rootlevel--;
+	}
+	return 0;
+    }
+    else {
+	return 1;
+    }
+}
+

lib/vector/rtree/indexm.c

@@ -0,0 +1,424 @@
+
+/****************************************************************************
+* 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 - file-based and memory-based R*-tree
+*               
+* PURPOSE:      Multidimensional index
+*
+* COPYRIGHT:    (C) 2010 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 <stdlib.h>
+#include <assert.h>
+#include <grass/gis.h>
+#include "index.h"
+#include "card.h"
+
+/* stack used for non-recursive insertion/deletion */
+struct stack
+{
+    struct Node *sn;	/* stack node */
+    int branch_id;	/* branch no to follow down */
+};
+
+int RTreeValidChildM(union Child *child)
+{
+    return (child->ptr != NULL);
+}
+
+/*
+ * Search in an index tree for all data retangles that
+ * overlap the argument rectangle.
+ * Return the number of qualifying data rects.
+ */
+int RTreeSearchM(struct RTree *t, struct Rect *r,
+                 SearchHitCallback *shcb, void *cbarg)
+{
+    struct Node *n;
+    int hitCount = 0, found;
+    int i;
+    struct stack s[MAXLEVEL];
+    int top = 0;
+
+    /* stack size of t->rootlevel + 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.ptr &&
+		    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.ptr;
+		    s[top].branch_id = 0;
+		    found = 0;
+		    break;
+		}
+	    }
+	    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.id &&
+		    RTreeOverlap(r, &(s[top].sn->branch[i].rect), t)) {
+		    hitCount++;
+		    if (shcb) {	/* call the user-provided callback */
+			if (!shcb(s[top].sn->branch[i].child.id, cbarg)) {
+			    /* callback wants to terminate search early */
+			    return hitCount;
+			}
+		    }
+		}
+	    }
+	    top--;
+	}
+    }
+
+    return hitCount;
+}
+
+/*
+ * Inserts a new data rectangle into the index structure.
+ * 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 RTreeInsertRect2M(struct Rect *r, union Child child, int level,
+			     struct Node **newnode,
+			     struct RTree *t,
+			     struct ListBranch **ee, int *overflow)
+{
+    int i;
+    struct Branch b;
+    struct Node *n, *n2;
+    struct Rect *cover;
+    struct stack s[MAXLEVEL];
+    int top = 0, down = 0;
+    int result;
+
+    /* 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.ptr;
+    }
+
+    /* Have reached level for insertion. Remove p rectangles or split */
+    if (s[top].sn->level == level) {
+	b.rect = *r;
+	/* child field of leaves contains tid of data record */
+	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);
+	/* update node count */
+	if (result == 1) {
+	    t->n_nodes++;
+	}
+    }
+    else {
+	/* Not supposed to happen */
+	assert(0);
+	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.ptr = n2;
+	    b.rect = RTreeNodeCover(b.child.ptr, 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);
+
+	    /* update node count */
+	    if (result == 1) {
+		t->n_nodes++;
+	    }
+	}
+    }
+
+    *newnode = n2;
+
+    return result;
+}
+
+/* 
+ * Insert a data rectangle into an index structure.
+ * RTreeInsertRectM 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 actual insertion.
+ */
+int RTreeInsertRectM(struct Rect *r, union Child child, int level,
+		     struct RTree *t)
+{
+    struct Node *newnode;
+    struct Node *newroot;
+    struct Branch b;
+    struct ListBranch *reInsertList = NULL;
+    struct ListBranch *e;
+    int result;
+    int i, overflow[MAXLEVEL];
+
+    /* R*-tree forced reinsertion: for each level only once */
+    for (i = 0; i < MAXLEVEL; i++)
+	overflow[i] = 1;
+
+    result =
+	RTreeInsertRect2M(r, child, level, &newnode, t, &reInsertList,
+			 overflow);
+			 
+    if (result == 1) {		/* root split */
+	/* grow a new root, & tree taller */
+	t->rootlevel++;
+	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);
+	/* branch for new node created by RTreeInsertRect2() */
+	b.rect = RTreeNodeCover(newnode, t);
+	b.child.ptr = newnode;
+	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 =
+		RTreeInsertRect2M(&(b.rect), b.child, level, &newnode, t,
+				 &reInsertList, overflow);
+
+	    if (result == 1) {	/* root split */
+		/* grow a new root, & tree taller */
+		t->rootlevel++;
+		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);
+		/* branch for new node created by RTreeInsertRect2() */
+		b.rect = RTreeNodeCover(newnode, t);
+		b.child.ptr = newnode; 
+		RTreeAddBranch(&b, newroot, NULL, NULL, NULL, NULL, t);
+		/* set new root node */
+		t->root = newroot;
+		t->n_nodes++;
+	    }
+	}
+    }
+
+    return result;
+}
+
+/*
+ * Delete a rectangle from non-root part of an index structure.
+ * 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
+RTreeDeleteRect2M(struct Rect *r, union Child child, struct RTree *t,
+		 struct ListNode **ee)
+{
+    int i, notfound = 1;
+    struct Node *n;
+    struct stack s[MAXLEVEL];
+    int top = 0, down = 0;
+    int minfill;
+
+    assert(ee);
+
+    /* 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.ptr &&
+		    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.ptr;
+		    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.id &&
+		    s[top].sn->branch[i].child.id == child.id) { /* found item */
+		    RTreeDisconnectBranch(s[top].sn, i, t);
+		    t->n_leafs--;
+		    notfound = 0;
+		    break;
+		}
+	    }
+	    if (notfound)	/* continue searching */
+		top--;
+	}
+    }
+
+    if (notfound) {
+	return notfound;
+    }
+
+    /* go back up */
+    while (top) {
+	down = top;
+	top--;
+	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.ptr, ee);
+	    RTreeDisconnectBranch(s[top].sn, i, t);
+	}
+    }
+
+    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 RTree.
+ * Returns 1 if record not found, 0 if success.
+ * RTreeDeleteRect1 provides for eliminating the root.
+ */
+int RTreeDeleteRectM(struct Rect *r, union Child child, struct RTree *t)
+{
+    int i;
+    struct Node *n;
+    struct ListNode *reInsertList = NULL;
+    struct ListNode *e;
+
+    if (!RTreeDeleteRect2M(r, child, t, &reInsertList)) {
+	/* found and deleted a data item */
+
+	/* reinsert any branches from eliminated nodes */
+	while (reInsertList) {
+	    t->n_nodes--;
+	    n = reInsertList->node;
+	    if (n->level > 0) {  /* reinsert node branches */
+		for (i = 0; i < t->nodecard; i++) {
+		    if (n->branch[i].child.ptr) {
+			RTreeInsertRectM(&(n->branch[i].rect),
+					 n->branch[i].child, n->level, t);
+		    }
+		}
+	    }
+	    else {  /* reinsert leaf branches */
+		for (i = 0; i < t->leafcard; i++) {
+		    if (n->branch[i].child.id) {
+			RTreeInsertRectM(&(n->branch[i].rect),
+					 n->branch[i].child, n->level, t);
+		    }
+		}
+	    }
+	    e = reInsertList;
+	    reInsertList = reInsertList->next;
+	    RTreeFreeNode(e->node);
+	    RTreeFreeListNode(e);
+	}
+
+	/* check for redundant root (not leaf, 1 child) and eliminate */
+	n = t->root;
+
+	if (n->count == 1 && n->level > 0) {
+	    for (i = 0; i < t->nodecard; i++) {
+		if (n->branch[i].child.ptr)
+		    break;
+	    }
+	    t->root = n->branch[i].child.ptr;
+	    RTreeFreeNode(n);
+	    t->rootlevel--;
+	}
+	return 0;
+    }
+    else {
+	return 1;
+    }
+}

lib/vector/rtree/io.c

@@ -0,0 +1,134 @@
+
+/****************************************************************************
+* 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 - file-based and memory-based R*-tree
+*               
+* 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 <stdlib.h>
+#include <stdio.h>
+#include <sys/types.h>
+#include <unistd.h>
+#include <assert.h>
+#include <grass/config.h>
+#include "index.h"
+
+/* add new free node position for recycling */
+void RTreeAddNodePos(off_t pos, int level, struct RTree *t)
+{
+    int which;
+    
+    if (t->free_nodes.avail >= t->free_nodes.alloc) {
+	size_t size;
+
+	t->free_nodes.alloc += 100;
+	size = t->free_nodes.alloc * sizeof(off_t);
+	t->free_nodes.pos = (off_t *)realloc((void *)t->free_nodes.pos, size);
+	assert(t->free_nodes.pos);
+    }
+    t->free_nodes.pos[t->free_nodes.avail++] = pos;
+    
+    which = pos == t->nb[level][1].pos;
+    t->nb[level][which].pos = -1;
+    t->nb[level][which].dirty = 0;
+    t->mru[level] = which == 0;
+}
+
+/* looks for free node position, sets file pointer, returns position */
+off_t RTreeGetNodePos(struct RTree *t)
+{
+    if (t->free_nodes.avail > 0) {
+	t->free_nodes.avail--;
+	return lseek(t->fd, t->free_nodes.pos[t->free_nodes.avail], SEEK_SET);
+    }
+    else {
+	return lseek(t->fd, 0, SEEK_END);
+    }
+}
+
+/* read node from file */
+size_t RTreeReadNode(struct Node *n, off_t nodepos, struct RTree *t)
+{
+    lseek(t->fd, nodepos, SEEK_SET);
+    return read(t->fd, n, t->nodesize);
+}
+
+/* get node from buffer or file */
+void RTreeGetNode(struct Node *n, off_t nodepos, int level, struct RTree *t)
+{
+    int which = nodepos == t->nb[level][1].pos;
+
+    if (t->nb[level][which].pos != nodepos) {
+	/* which is 0 */
+	/* replace least recently used */
+	/* least recently used is faster than most recently used */
+	if (t->nb[level][which].pos != -1)
+	    which = t->mru[level] == 0;
+	/* rewrite node in buffer */
+	if (t->nb[level][which].dirty) {
+	    RTreeRewriteNode(&(t->nb[level][which].n), t->nb[level][which].pos, t);
+	    t->nb[level][which].dirty = 0;
+	}
+	RTreeReadNode(&(t->nb[level][which].n), nodepos, t);
+	t->nb[level][which].pos = nodepos;
+    }
+    t->mru[level] = which;
+    *n = t->nb[level][which].n;
+}
+
+/* write new node to file */
+size_t RTreeWriteNode(struct Node *n, struct RTree *t)
+{
+    /* file position must be set first with RTreeGetFNodePos() */
+    return write(t->fd, n, t->nodesize);
+}
+
+/* rewrite updated node to file */
+size_t RTreeRewriteNode(struct Node *n, off_t nodepos, struct RTree *t)
+{
+    lseek(t->fd, nodepos, SEEK_SET);
+    return write(t->fd, n, t->nodesize);
+}
+
+/* update node in buffer */
+void RTreePutNode(struct Node *n, off_t nodepos, struct RTree *t)
+{
+    int which = nodepos == t->nb[n->level][1].pos;
+    
+    t->nb[n->level][which].n = *n;
+    t->nb[n->level][which].dirty = 1;
+    t->mru[n->level] = which;
+}
+
+/* update rectangle */
+void RTreeUpdateRect(struct Rect *r, struct Node *n, off_t nodepos, int b, struct RTree *t)
+{
+    int which = nodepos == t->nb[n->level][1].pos;
+    
+    t->nb[n->level][which].n.branch[b].rect = n->branch[b].rect = *r;
+    t->nb[n->level][which].dirty = 1;
+    t->mru[n->level] = which;
+}
+
+void RTreeFlushBuffer(struct RTree *t)
+{
+    int i;
+    
+    for (i = 0; i <= t->rootlevel; i++) {
+	if (t->nb[i][0].dirty)
+	    RTreeRewriteNode(&(t->nb[i][0].n), t->nb[i][0].pos, t);
+	if (t->nb[i][1].dirty)
+	    RTreeRewriteNode(&(t->nb[i][1].n), t->nb[i][1].pos, t);
+    }
+}

lib/vector/rtree/node.c

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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
@@ -19,6 +19,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <assert.h>
+#include <grass/gis.h>
 #include "index.h"
 #include "card.h"
 
@@ -30,11 +31,17 @@ struct dist
 };
 
 /* Initialize one branch cell in an internal node. */
-static void RTreeInitNodeBranch(struct Branch *b)
+static void RTreeInitNodeBranchM(struct Branch *b)
 {
     RTreeInitRect(&(b->rect));
     b->child.ptr = NULL;
+}
 
+/* Initialize one branch cell in an internal node. */
+static void RTreeInitNodeBranchF(struct Branch *b)
+{
+    RTreeInitRect(&(b->rect));
+    b->child.pos = -1;
 }
 
 /* Initialize one branch cell in a leaf node. */
@@ -42,28 +49,23 @@ static void RTreeInitLeafBranch(struct Branch *b)
 {
     RTreeInitRect(&(b->rect));
     b->child.id = 0;
-
 }
 
-static void (*RTreeInitBranch[2]) () = {
-    RTreeInitLeafBranch, RTreeInitNodeBranch
+static void (*RTreeInitBranch[3]) () = {
+    RTreeInitLeafBranch, RTreeInitNodeBranchM, RTreeInitNodeBranchF
 };
 
 /* Initialize a Node structure. */
-void RTreeInitNode(struct Node *n, int level)
+/* type = 1: leaf, type = 2: internal, memory, type = 3: internal, file */
+void RTreeInitNode(struct Node *n, int type)
 {
     int i;
 
     n->count = 0;
     n->level = -1;
-    if (level > 0) {
-	for (i = 0; i < MAXCARD; i++)
-	    RTreeInitNodeBranch(&(n->branch[i]));
-    }
-    else {
-	for (i = 0; i < MAXCARD; i++)
-	    RTreeInitLeafBranch(&(n->branch[i]));
-    }
+
+    for (i = 0; i < MAXCARD; i++)
+	RTreeInitBranch[type](&(n->branch[i]));
 }
 
 /* Make a new node and initialize to have all branch cells empty. */
@@ -73,7 +75,7 @@ struct Node *RTreeNewNode(struct RTree *t, int level)
 
     n = (struct Node *)malloc((size_t) t->nodesize);
     assert(n);
-    RTreeInitNode(n, level);
+    RTreeInitNode(n, NODETYPE(level, t->fd));
     return n;
 }
 
@@ -92,12 +94,10 @@ struct Rect RTreeNodeCover(struct Node *n, struct RTree *t)
     int i, first_time = 1;
     struct Rect r;
 
-    assert(n);
-
     RTreeInitRect(&r);
     if ((n)->level > 0) { /* internal node */
 	for (i = 0; i < t->nodecard; i++) {
-	    if (n->branch[i].child.ptr) {
+	    if (t->valid_child(&(n->branch[i].child))) {
 		if (first_time) {
 		    r = n->branch[i].rect;
 		    first_time = 0;
@@ -134,23 +134,21 @@ struct Rect RTreeNodeCover(struct Node *n, struct RTree *t)
  * to get the best resolution when searching.
  */
 
-static int RTreePickBranch1(struct Rect *r, struct Node *n, struct RTree *t)
+static int RTreePickLeafBranch(struct Rect *r, struct Node *n, struct RTree *t)
 {
     struct Rect *rr;
     int i, j;
-    RectReal increase, bestIncr = (RectReal) - 1, area, bestArea = 0;
+    RectReal increase, bestIncr = -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.ptr) {
+	if (t->valid_child(&(n->branch[i].child))) {
 	    rr = &n->branch[i].rect;
 	    tmp_rect = RTreeCombineRect(r, rr, t);
 	    area = RTreeRectSphericalVolume(rr, t);
@@ -199,18 +197,17 @@ int RTreePickBranch(struct Rect *r, struct Node *n, struct RTree *t)
 {
     struct Rect *rr;
     int i, first_time = 1;
-    RectReal increase, bestIncr = (RectReal) - 1, area, bestArea = 0;
+    RectReal increase, bestIncr = (RectReal) -1, area, bestArea = 0;
     int best = 0;
     struct Rect tmp_rect;
 
-    assert(r && n && t);
     assert((n)->level > 0);	/* must not be called on leaf node */
 
     if ((n)->level == 1)
-	return RTreePickBranch1(r, n, t);
+	return RTreePickLeafBranch(r, n, t);
 
     for (i = 0; i < t->nodecard; i++) {
-	if (n->branch[i].child.ptr) {
+	if (t->valid_child(&(n->branch[i].child))) {
 	    rr = &n->branch[i].rect;
 	    area = RTreeRectSphericalVolume(rr, t);
 	    tmp_rect = RTreeCombineRect(r, rr, t);
@@ -235,8 +232,11 @@ void RTreeDisconnectBranch(struct Node *n, int i, struct RTree *t)
 {
     if ((n)->level > 0) {
 	assert(n && i >= 0 && i < t->nodecard);
-	assert(n->branch[i].child.ptr);
-	RTreeInitNodeBranch(&(n->branch[i]));
+	assert(t->valid_child(&(n->branch[i].child)));
+	if (t->fd < 0)
+	    RTreeInitNodeBranchM(&(n->branch[i]));
+	else
+	    RTreeInitNodeBranchF(&(n->branch[i]));
     }
     else {
 	assert(n && i >= 0 && i < t->leafcard);
@@ -265,11 +265,11 @@ void RTreeDestroyNode(struct Node *n, int nodes)
     RTreeFreeNode(n);
 }
 
-/**********************************************************************
- *                                                                    *
- *   R*-tree: force remove p (currently 3) branches for reinsertion   *
- *                                                                    *
- **********************************************************************/
+/****************************************************************
+ *                                                              *
+ *   R*-tree: force remove FORCECARD branches for reinsertion   *
+ *                                                              *
+ ****************************************************************/
 
 /*
  * swap dist structs
@@ -399,15 +399,8 @@ static void RTreeQuicksortDist(struct dist *d, int n)
 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.
@@ -424,11 +417,16 @@ static void RTreeReInsertBranch(struct Branch b, int level,
     *ee = l;
 }
 
+/*
+ * Remove branches from a node. Select the 2 branches whose rectangle 
+ * center is farthest away from node cover center.
+ * Old node updated.
+ */
 static void RTreeRemoveBranches(struct Node *n, struct Branch *b,
 				struct ListBranch **ee, struct Rect *cover,
 				struct RTree *t)
 {
-    int i, j, maxkids, is_node;
+    int i, j, maxkids, type;
     RectReal center_n[NUMDIMS], center_r, delta;
     struct Branch branchbuf[MAXCARD + 1];
     struct dist rdist[MAXCARD + 1];
@@ -436,15 +434,9 @@ static void RTreeRemoveBranches(struct Node *n, struct Branch *b,
 
     assert(cover);
 
-    if ((n)->level > 0) {
-	maxkids = t->nodecard;
-	is_node = 1;
-    }
-    else {
-	maxkids = t->leafcard;
-	is_node = 0;
-    }
-
+    maxkids = MAXKIDS((n)->level, t);
+    type = NODETYPE((n)->level, t->fd);
+    
     assert(n->count == maxkids);	/* must be full */
 
     new_cover = RTreeCombineRect(cover, &(b->rect), t);
@@ -467,7 +459,7 @@ static void RTreeRemoveBranches(struct Node *n, struct Branch *b,
 	    rdist[i].distance += delta * delta;
 	}
 
-	RTreeInitBranch[is_node](&(n->branch[i]));
+	RTreeInitBranch[type](&(n->branch[i]));
     }
 
     /* new branch */
@@ -501,23 +493,20 @@ static void RTreeRemoveBranches(struct Node *n, struct Branch *b,
  * 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.
- * Returns 2 if branches wereremoved for forced reinsertion
+ * Returns 2 if branches were removed for forced reinsertion
  */
 int RTreeAddBranch(struct Branch *b, struct Node *n,
-		   struct Node **new_node, struct ListBranch **ee,
+		   struct Node **newnode, struct ListBranch **ee,
 		   struct Rect *cover, int *overflow, struct RTree *t)
 {
     int i, maxkids;
 
-    assert(b);
-    assert(n);
-
-    maxkids = ((n)->level > 0 ? t->nodecard : t->leafcard);
+    maxkids = MAXKIDS((n)->level, t);
 
     if (n->count < maxkids) {	/* split won't be necessary */
 	if ((n)->level > 0) {   /* internal node */
-	    for (i = 0; i < t->nodecard; i++) {	/* find empty branch */
-		if (n->branch[i].child.ptr == NULL) {
+	    for (i = 0; i < maxkids; i++) {	/* find empty branch */
+		if (!t->valid_child(&(n->branch[i].child))) {
 		    n->branch[i] = *b;
 		    n->count++;
 		    break;
@@ -526,7 +515,7 @@ int RTreeAddBranch(struct Branch *b, struct Node *n,
 	    return 0;
 	}
 	else if ((n)->level == 0) {   /* leaf */
-	    for (i = 0; i < t->leafcard; i++) {	/* find empty branch */
+	    for (i = 0; i < maxkids; i++) {	/* find empty branch */
 		if (n->branch[i].child.id == 0) {
 		    n->branch[i] = *b;
 		    n->count++;
@@ -537,15 +526,18 @@ int RTreeAddBranch(struct Branch *b, struct Node *n,
 	}
     }
     else {
-	if (n->level < t->n_levels && overflow[n->level]) {
+	if (n->level < t->rootlevel && overflow[n->level]) {
 	    /* R*-tree forced reinsert */
 	    RTreeRemoveBranches(n, b, ee, cover, t);
 	    overflow[n->level] = 0;
 	    return 2;
 	}
 	else {
-	    *new_node = RTreeNewNode(t, (n)->level);
-	    RTreeSplitNode(n, b, *new_node, t);
+	    if (t->fd > -1)
+		RTreeInitNode(*newnode, NODETYPE(n->level, t->fd));
+	    else
+		*newnode = RTreeNewNode(t, (n)->level);
+	    RTreeSplitNode(n, b, *newnode, t);
 	    return 1;
 	}
     }
@@ -580,7 +572,6 @@ void RTreePrintNode(struct Node *n, int depth, struct RTree *t)
 
     RTreeTabIn(depth);
 
-
     maxkids = (n->level > 0 ? t->nodecard : t->leafcard);
 
     fprintf(stdout, "node");
@@ -606,3 +597,4 @@ void RTreePrintNode(struct Node *n, int depth, struct RTree *t)
 	}
     }
 }
+

lib/vector/rtree/rect.c

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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
@@ -251,7 +251,7 @@ RectReal RTreeRectSphericalVolume(struct Rect *R, struct RTree *t)
 /*
  * The exact volume of the bounding sphere for the given Rect.
  */
-RectReal RTreeRectSphericalVolume(struct Rect * r, struct RTree * t)
+RectReal RTreeRectSphericalVolume(struct Rect *r, struct RTree *t)
 {
     int i;
     double sum_of_squares = 0, radius, half_extent;
@@ -272,7 +272,7 @@ RectReal RTreeRectSphericalVolume(struct Rect * r, struct RTree * t)
 /*-----------------------------------------------------------------------------
 | Calculate the n-dimensional surface area of a rectangle
 -----------------------------------------------------------------------------*/
-RectReal RTreeRectSurfaceArea(struct Rect * r, struct RTree * t)
+RectReal RTreeRectSurfaceArea(struct Rect *r, struct RTree *t)
 {
     int i, j;
     RectReal j_extent, sum = (RectReal) 0;
@@ -301,7 +301,7 @@ RectReal RTreeRectSurfaceArea(struct Rect * r, struct RTree * t)
 | 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)
+RectReal RTreeRectMargin(struct Rect *r, struct RTree *t)
 {
     int i;
     RectReal margin = 0.0;
@@ -342,6 +342,26 @@ struct Rect RTreeCombineRect(struct Rect *r, struct Rect *rr, struct RTree *t)
 
 
 /*-----------------------------------------------------------------------------
+| Decide whether two rectangles are identical.
+-----------------------------------------------------------------------------*/
+int RTreeCompareRect(struct Rect *r, struct Rect *s, struct RTree *t)
+{
+    register int i, j;
+
+    assert(r && s);
+
+    for (i = 0; i < t->ndims; i++) {
+	j = i + NUMDIMS;	/* index for high sides */
+	if (r->boundary[i] != s->boundary[i] ||
+	    r->boundary[j] != s->boundary[j]) {
+	    return 0;
+	}
+    }
+    return 1;
+}
+
+
+/*-----------------------------------------------------------------------------
 | Decide whether two rectangles overlap.
 -----------------------------------------------------------------------------*/
 int RTreeOverlap(struct Rect *r, struct Rect *s, struct RTree *t)
@@ -368,7 +388,7 @@ int RTreeContained(struct Rect *r, struct Rect *s, struct RTree *t)
 {
     register int i, j, result;
 
-    assert(r && s);		/* same as in RTreeOverlap() */
+    assert(r && s);
 
     /* undefined rect is contained in any other */
     if (Undefined(r))

lib/vector/rtree/rtree.h

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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/split.c

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* 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 
@@ -34,7 +34,7 @@ struct Rect CoverSplit;
 RectReal CoverSplitArea;
 
 /* variables for finding a partition */
-struct PartitionVars Partitions[METHODS];
+struct PartitionVars Partitions[1];
 
 /*----------------------------------------------------------------------
 | Load branch buffer with branches from full node plus the extra branch.
@@ -44,14 +44,11 @@ static void RTreeGetBranches(struct Node *n, struct Branch *b,
 {
     int i, maxkids = 0;
 
-    assert(n);
-    assert(b);
-
     if ((n)->level > 0) {
 	maxkids = t->nodecard;
 	/* load the branch buffer */
 	for (i = 0; i < maxkids; i++) {
-	    assert(n->branch[i].child.ptr);	/* n should have every entry full */
+	    assert(t->valid_child(&(n->branch[i].child)));	/* n should have every entry full */
 	    BranchBuf[i] = n->branch[i];
 	}
     }
@@ -67,14 +64,16 @@ static void RTreeGetBranches(struct Node *n, struct Branch *b,
     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);
+    if (METHOD == 0) { /* quadratic split */
+	/* 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);
     }
-    CoverSplitArea = RTreeRectSphericalVolume(&CoverSplit, t);
 
-    RTreeInitNode(n, n->level);
+    RTreeInitNode(n, NODETYPE(n->level, t->fd));
 }
 
 /*----------------------------------------------------------------------
@@ -83,26 +82,27 @@ static void RTreeGetBranches(struct Node *n, struct Branch *b,
 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);
+    if (METHOD == 0) {
+	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 Guttman's quadratic splitting method              *
- *                                                                    *
- **********************************************************************/
+/***************************************************
+ *                                                 *
+ *    Toni Guttman's quadratic splitting method    *
+ *                                                 *
+ ***************************************************/
 
 /*----------------------------------------------------------------------
 | Pick two rects from set to be the first elements of the two groups.
@@ -146,11 +146,6 @@ static void RTreeLoadNodes(struct Node *n, struct Node *q,
 {
     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)
@@ -167,8 +162,6 @@ 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;
@@ -188,8 +181,6 @@ 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);
@@ -242,8 +233,6 @@ static void RTreeMethodZero(struct PartitionVars *p, int minfill,
     RectReal biggestDiff;
     int group, chosen = 0, betterGroup = 0;
 
-    assert(p);
-
     RTreeInitPVars(p, BranchCount, minfill);
     RTreePickSeeds(p, t);
 
@@ -367,9 +356,9 @@ static int RTreePartitionBranchBuf(int first, int last, int side)
 	return last;
     }
 
-    /* Larger of two */
-    if (RTreeCompareBranches(&(BranchBuf[first]), &(BranchBuf[mid]), side) ==
-	1) {
+    /* larger of two */
+    if (RTreeCompareBranches(&(BranchBuf[first]), &(BranchBuf[mid]), side)
+	== 1) {
 	larger = pivot = first;
 	smaller = mid;
     }
@@ -462,9 +451,9 @@ static void RTreeQuicksortBranchBuf(int side)
 | 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 maxkids, struct RTree *t)
 {
-    int i, j, k, l, s, maxkids;
+    int i, j, k, l, s;
     int axis = 0, best_axis = 0, side = 0, best_side[NUMDIMS];
     int best_cut[NUMDIMS];
     RectReal margin, smallest_margin = 0;
@@ -472,12 +461,8 @@ static void RTreeMethodOne(struct PartitionVars *p, int minfill,
     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);
-
     margin = DBL_MAX;
 
     /* choose split axis */
@@ -523,8 +508,8 @@ static void RTreeMethodOne(struct PartitionVars *p, int minfill,
 
 		/* the margin is the sum of the lengths of the edges of a rectangle */
 		margin =
-		    RTreeRectMargin(&testrect1,
-				    t) + RTreeRectMargin(&testrect2, t);
+		    RTreeRectMargin(&testrect1, t) +
+		    RTreeRectMargin(&testrect2, t);
 
 		/* remember best axis */
 		if (margin <= smallest_margin) {
@@ -539,10 +524,8 @@ static void RTreeMethodOne(struct 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 (testrect1.boundary[k] > testrect2.boundary[k + NUMDIMS] ||
+			testrect1.boundary[k + NUMDIMS] < testrect2.boundary[k]) {
 			overlap = 0;
 			break;
 		    }
@@ -564,8 +547,8 @@ static void RTreeMethodOne(struct PartitionVars *p, int minfill,
 		    overlap = RTreeRectVolume(&orect, t);
 
 		vol =
-		    RTreeRectVolume(&testrect1,
-				    t) + RTreeRectVolume(&testrect2, t);
+		    RTreeRectVolume(&testrect1, t) +
+		    RTreeRectVolume(&testrect2, t);
 
 		/* get best cut for this axis */
 		if (overlap <= smallest_overlap) {
@@ -582,9 +565,9 @@ static void RTreeMethodOne(struct PartitionVars *p, int minfill,
 			best_side[i] = s;
 		    }
 		}
-	    }			/* end of distribution check */
-	}			/* end of side check */
-    }				/* end of axis check */
+	    }    /* 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)
@@ -614,18 +597,17 @@ 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, 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);
+    
+    if (METHOD == 1) /* R* split */
+	RTreeMethodOne(p, MINFILL(level, t), MAXKIDS(level, t), t);
+    else
+	RTreeMethodZero(p, MINFILL(level, t), t);
 
     /*
      * put branches from buffer into 2 nodes

lib/vector/rtree/split.h

@@ -5,11 +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
+*               Markus Metz - file-based and memory-based R*-tree
 *               
 * PURPOSE:      Multidimensional index
 *
-* COPYRIGHT:    (C) 2009 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 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
@@ -20,7 +20,9 @@
 | Definitions and global variables.
 -----------------------------------------------------------------------------*/
 
-#define METHODS 1
+/* METHOD 0: R-Tree, quadratic split */
+/* METHOD 1: R*-Tree split */
+#define METHOD 1
 
 struct PartitionVars {
     int partition[MAXCARD + 1];
@@ -37,5 +39,6 @@ extern struct Rect CoverSplit;
 extern RectReal CoverSplitArea;
 
 /* variables for finding a partition */
-extern struct PartitionVars Partitions[METHODS];
+extern struct PartitionVars Partitions[1];
 
+extern void RTreeInitPVars(struct PartitionVars *, int, int);