btree2 lib: improved k-d tree balancing

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

lib/btree2/kdtree.c

@@ -36,7 +36,7 @@ static int rcallsmax = 0;
 
 static struct kdnode *kdtree_insert2(struct kdtree *, struct kdnode *,
                                      struct kdnode *, int, int);
-static int kdtree_replace(struct kdtree *, struct kdnode *, int);
+static int kdtree_replace(struct kdtree *, struct kdnode *);
 static int kdtree_balance(struct kdtree *, struct kdnode *, int);
 static int kdtree_first(struct kdtrav *, double *, int *);
 static int kdtree_next(struct kdtrav *, double *, int *);
@@ -70,6 +70,7 @@ static struct kdnode *kdtree_newnode(struct kdtree *t)
     n->c = G_malloc(t->ndims * sizeof(double));
     n->dim = 0;
     n->depth = 0;
+    n->balance = 0;
     n->uid = 0;
     n->child[0] = NULL;
     n->child[1] = NULL;
@@ -83,6 +84,32 @@ static void kdtree_free_node(struct kdnode *n)
     G_free(n);
 }
 
+static void kdtree_update_node(struct kdtree *t, struct kdnode *n)
+{
+    int ld, rd, btol;
+
+    ld = (!n->child[0] ? -1 : n->child[0]->depth);
+    rd = (!n->child[1] ? -1 : n->child[1]->depth);
+    n->depth = MAX(ld, rd) + 1;
+
+    n->balance = 0;
+    /* set balance flag if any of the node's subtrees needs balancing
+     * or if the node itself needs balancing */
+    if ((n->child[0] && n->child[0]->balance) ||
+        (n->child[1] && n->child[1]->balance)) {
+	n->balance = 1;
+
+	return;
+    }
+
+    btol = t->btol;
+    if (!n->child[0] || !n->child[1])
+	btol = 2;
+
+    if (ld > rd + btol || rd > ld + btol) 
+	n->balance = 1;
+}
+
 /* create a new k-d tree with ndims dimensions,
  * optionally set balancing tolerance */
 struct kdtree *kdtree_create(char ndims, int *btol)
@@ -185,7 +212,7 @@ int kdtree_remove(struct kdtree *t, double *c, int uid)
     } s[256];
     int top;
     int dir, found;
-    int ld, rd;
+    int balance, bmode;
 
     sn.c = c;
     sn.uid = uid;
@@ -221,8 +248,8 @@ int kdtree_remove(struct kdtree *t, double *c, int uid)
 	    dir = s[top].dir;
 	    n->child[dir] = NULL;
 
-	    /* update node depth */
-	    n->depth = (!n->child[!dir] ? 0 : n->child[!dir]->depth + 1);
+	    /* update node */
+	    kdtree_update_node(t, n);
 	}
 	else {
 	    t->root = NULL;
@@ -231,30 +258,73 @@ int kdtree_remove(struct kdtree *t, double *c, int uid)
 	}
     }
     else
-	kdtree_replace(t, s[top].n, 1);
+	kdtree_replace(t, s[top].n);
 
-    if (top) {
+    while (top) {
 	top--;
-	dir = s[top].dir;
 	n = s[top].n;
 
-	while (kdtree_balance(t, n->child[dir], 0));
-
-	/* update node depth */
-	ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	n->depth = MAX(ld, rd) + 1;
+	/* update node */
+	kdtree_update_node(t, n);
     }
-    while (top) {
-	top--;
-	n = s[top].n;
 
-	/* update node depth */
-	ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	n->depth = MAX(ld, rd) + 1;
+    balance = 1;
+    bmode = 1;
+    if (balance) {
+	struct kdnode *r;
+	int iter, bmode2;
 
-	while (kdtree_balance(t, n, 0));
+	/* fix any inconsistencies in the (sub-)tree */
+	iter = 0;
+	bmode2 = 0;
+	top = 0;
+	r = t->root;
+	s[top].n = r;
+	while (top >= 0) {
+
+	    n = s[top].n;
+
+	    /* top-down balancing
+	     * slower but more compact */
+	    if (!bmode2) {
+		while (kdtree_balance(t, n, bmode));
+	    }
+
+	    /* go down */
+	    if (n->child[0] && n->child[0]->balance) {
+		dir = 0;
+		top++;
+		s[top].n = n->child[dir];
+	    }
+	    else if (n->child[1] && n->child[1]->balance) {
+		dir = 1;
+		top++;
+		s[top].n = n->child[dir];
+	    }
+	    /* go back up */
+	    else {
+
+		/* bottom-up balancing
+		 * faster but less compact */
+		kdtree_update_node(t, n);
+		if (bmode2) {
+		    while (kdtree_balance(t, n, bmode));
+		}
+		top--;
+		if (top >= 0) {
+		    kdtree_update_node(t, s[top].n);
+		}
+		if (!bmode2 && top == 0) {
+		    iter++;
+		    if (iter == 2) {
+			/* the top node has been visited twice, 
+			 * switch from top-down to bottom-up balancing */
+			iter = 0;
+			bmode2 = 1;
+		    }
+		}
+	    }
+	}
     }
 
     return 1;
@@ -291,6 +361,15 @@ void kdtree_optimize(struct kdtree *t, int level)
 
 	ld = (!n->child[0] ? -1 : n->child[0]->depth);
 	rd = (!n->child[1] ? -1 : n->child[1]->depth);
+	
+	if (ld < rd)
+	    while (kdtree_balance(t, n->child[0], level));
+	else if (ld > rd)
+	    while (kdtree_balance(t, n->child[1], level));
+
+	ld = (!n->child[0] ? -1 : n->child[0]->depth);
+	rd = (!n->child[1] ? -1 : n->child[1]->depth);
+	n->depth = MAX(ld, rd) + 1;
 
 	dir = (rd > ld);
 
@@ -303,18 +382,17 @@ void kdtree_optimize(struct kdtree *t, int level)
 	n = s[top].n;
 
 	/* balance node */
+	while (kdtree_balance(t, n, level)) {
+	    nbal++;
+	}
 	while (kdtree_balance(t, n->child[0], level));
 	while (kdtree_balance(t, n->child[1], level));
+
 	ld = (!n->child[0] ? -1 : n->child[0]->depth);
 	rd = (!n->child[1] ? -1 : n->child[1]->depth);
 	n->depth = MAX(ld, rd) + 1;
-	while (kdtree_balance(t, n, level)) {
-	    while (kdtree_balance(t, n->child[0], level));
-	    while (kdtree_balance(t, n->child[1], level));
 
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
+	while (kdtree_balance(t, n, level)) {
 	    nbal++;
 	}
     }
@@ -323,18 +401,17 @@ void kdtree_optimize(struct kdtree *t, int level)
 	n = s[top].n;
 
 	/* balance node */
+	while (kdtree_balance(t, n, level)) {
+	    nbal++;
+	}
 	while (kdtree_balance(t, n->child[0], level));
 	while (kdtree_balance(t, n->child[1], level));
+
 	ld = (!n->child[0] ? -1 : n->child[0]->depth);
 	rd = (!n->child[1] ? -1 : n->child[1]->depth);
 	n->depth = MAX(ld, rd) + 1;
-	while (kdtree_balance(t, n, level)) {
-	    while (kdtree_balance(t, n->child[0], level));
-	    while (kdtree_balance(t, n->child[1], level));
 
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
+	while (kdtree_balance(t, n, level)) {
 	    nbal++;
 	}
 
@@ -364,19 +441,17 @@ void kdtree_optimize(struct kdtree *t, int level)
 	    n = s[top].n;
 
 	    /* balance node */
+	    while (kdtree_balance(t, n, level)) {
+		nbal++;
+	    }
 	    while (kdtree_balance(t, n->child[0], level));
 	    while (kdtree_balance(t, n->child[1], level));
 
 	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
 	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
 	    n->depth = MAX(ld, rd) + 1;
-	    while (kdtree_balance(t, n, level)) {
-		while (kdtree_balance(t, n->child[0], level));
-		while (kdtree_balance(t, n->child[1], level));
 
-		ld = (!n->child[0] ? -1 : n->child[0]->depth);
-		rd = (!n->child[1] ? -1 : n->child[1]->depth);
-		n->depth = MAX(ld, rd) + 1;
+	    while (kdtree_balance(t, n, level)) {
 		nbal++;
 	    }
 
@@ -417,7 +492,7 @@ void kdtree_optimize(struct kdtree *t, int level)
 	}
     }
 
-    G_debug(1, "k-d tree optimization: %d times balanced, new depth %d",
+    G_debug(0, "k-d tree optimization: %d times balanced, new depth %d",
             nbal, t->root->depth);
 
     return;
@@ -801,11 +876,11 @@ int kdtree_traverse(struct kdtrav *trav, double *c, int *uid)
 /*            internal functions              */
 /**********************************************/
 
-static int kdtree_replace(struct kdtree *t, struct kdnode *r, int bmode)
+static int kdtree_replace(struct kdtree *t, struct kdnode *r)
 {
     double mindist;
     int rdir, ordir, dir;
-    int ld, rd, old_depth;
+    int ld, rd;
     struct kdnode *n, *rn, *or;
     struct kdstack {
 	struct kdnode *n;
@@ -1013,24 +1088,8 @@ static int kdtree_replace(struct kdtree *t, struct kdnode *r, int bmode)
     n->child[dir] = NULL;
     t->count--;
 
-    old_depth = n->depth;
-
-    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-    n->depth = MAX(ld, rd) + 1;
-
-#if 0
-    /* this can cause stack overflow because of too high levels of recursion */
-    if (bmode > 1)
-	while (kdtree_balance(t, n, bmode));
-#endif
-
-    if (n->depth == old_depth)
-	top = 0;
-
-#ifdef KD_DEBUG
-    top = top2 - 1;
-#endif
+    kdtree_update_node(t, n);
+    top++;
 
     /* go back up */
     while (top) {
@@ -1049,10 +1108,8 @@ static int kdtree_replace(struct kdtree *t, struct kdnode *r, int bmode)
 	}
 #endif
 
-	/* update depth */
-	ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	n->depth = MAX(ld, rd) + 1;
+	/* update node */
+	kdtree_update_node(t, n);
     }
 
     return nr;
@@ -1076,7 +1133,7 @@ static int kdtree_balance(struct kdtree *t, struct kdnode *r, int bmode)
     
     if (old_depth != r->depth) {
 	G_warning("balancing: depth is wrong: %d != %d", r->depth, old_depth);
-	r->depth = old_depth;
+	kdtree_update_node(t, r);
     }
 
     /* subtree difference */
@@ -1101,7 +1158,7 @@ static int kdtree_balance(struct kdtree *t, struct kdnode *r, int bmode)
     or->uid = r->uid;
     or->dim = t->nextdim[r->dim];
 
-    if (!kdtree_replace(t, r, bmode))
+    if (!kdtree_replace(t, r))
 	G_fatal_error("kdtree_balance: nothing replaced");
 
 #ifdef KD_DEBUG
@@ -1113,12 +1170,10 @@ static int kdtree_balance(struct kdtree *t, struct kdnode *r, int bmode)
     }
 #endif
 
-    r->child[!dir] = kdtree_insert2(t, r->child[!dir], or, bmode, 1);
+    r->child[!dir] = kdtree_insert2(t, r->child[!dir], or, bmode, 1); /* bmode */
 
-    /* update node depth */
-    ld = (!r->child[0] ? -1 : r->child[0]->depth);
-    rd = (!r->child[1] ? -1 : r->child[1]->depth);
-    r->depth = MAX(ld, rd) + 1;
+    /* update node */
+    kdtree_update_node(t, r);
 
     if (r->depth == old_depth) {
 	G_debug(4, "balancing had no effect");
@@ -1142,9 +1197,6 @@ static struct kdnode *kdtree_insert2(struct kdtree *t, struct kdnode *r,
     } s[256];
     int top;
     int dir;
-    int ld, rd;
-    int old_depth;
-    int go_back;
     int bmode;
 
     if (!r) {
@@ -1160,46 +1212,20 @@ static struct kdnode *kdtree_insert2(struct kdtree *t, struct kdnode *r,
 	rcallsmax = rcalls;
 
     /* balancing modes
-     * bmode = 0: a posteriori, top-down and bottom-up balancing, no recursion
-     *            slow, high tree depth (sub-optimal tree)
-     * bmode = 1: top-down and bottom-up balancing, recursion
-     *            faster, lower tree depth (more optimal tree)
-     * bmode = 2: bottom-up balancing
-     *            fastest, tree depth similar to bmode = 1
-     * most optimal, but slowest: recursion, bottom-up, a posteriori
+     * bmode = 0: no recursion (only insert -> balance -> insert)
+     *            slower, higher tree depth
+     * bmode = 1: recursion (insert -> balance -> insert -> balance ...)
+     *            faster, more compact tree
      *  */
-    bmode = 2;
+    bmode = 1;
 
     /* find node with free child */
     top = 0;
-    go_back = 0;
     s[top].n = r;
     while (s[top].n) {
 
 	n = s[top].n;
 
-	if (balance && bmode < 2) {
-	    /* balance node */
-	    old_depth = n->depth;
-
-	    while (kdtree_balance(t, n->child[0], bmode));
-	    while (kdtree_balance(t, n->child[1], bmode));
-
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
-	    while (kdtree_balance(t, n, bmode)) {
-		while (kdtree_balance(t, n->child[0], bmode));
-		while (kdtree_balance(t, n->child[1], bmode));
-
-		ld = (!n->child[0] ? -1 : n->child[0]->depth);
-		rd = (!n->child[1] ? -1 : n->child[1]->depth);
-		n->depth = MAX(ld, rd) + 1;
-	    }
-	    if (old_depth != n->depth)
-		go_back = top;
-	}
-
 	if (!cmpc(nnew, n, t) && (!dc || nnew->uid == n->uid)) {
 
 	    G_debug(1, "KD node exists already, nothing to do");
@@ -1230,30 +1256,7 @@ static struct kdnode *kdtree_insert2(struct kdtree *t, struct kdnode *r,
 	nnew->dim = t->nextdim[n->dim];
 
 	t->count++;
-
-	old_depth = n->depth;
-	n->depth = (!n->child[!dir] ? 1 : n->child[!dir]->depth + 1);
-
-	if (balance) {
-	    /* balance parent */
-	    while (kdtree_balance(t, n->child[0], bmode));
-	    while (kdtree_balance(t, n->child[1], bmode));
-
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
-	    while (kdtree_balance(t, n, bmode)) {
-		while (kdtree_balance(t, n->child[0], bmode));
-		while (kdtree_balance(t, n->child[1], bmode));
-
-		ld = (!n->child[0] ? -1 : n->child[0]->depth);
-		rd = (!n->child[1] ? -1 : n->child[1]->depth);
-		n->depth = MAX(ld, rd) + 1;
-	    }
-	}
-
-	if (old_depth != n->depth)
-	    go_back = top;
+	top++;
     }
 
     /* go back up */
@@ -1261,31 +1264,10 @@ static struct kdnode *kdtree_insert2(struct kdtree *t, struct kdnode *r,
 	top--;
 	n = s[top].n;
 
-	/* update node depth */
-	ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	n->depth = MAX(ld, rd) + 1;
-
-	if (balance) {
-	    /* balance node */
-	    /* slightly reduced tree depth, a bit slower: */
-	    /*
-	    while (kdtree_balance(t, n->child[0], bmode));
-	    while (kdtree_balance(t, n->child[1], bmode));
-	    */
+	/* update node */
+	kdtree_update_node(t, n);
 
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
-	    while (kdtree_balance(t, n, bmode)) {
-		while (kdtree_balance(t, n->child[0], bmode));
-		while (kdtree_balance(t, n->child[1], bmode));
-
-		ld = (!n->child[0] ? -1 : n->child[0]->depth);
-		rd = (!n->child[1] ? -1 : n->child[1]->depth);
-		n->depth = MAX(ld, rd) + 1;
-	    }
-	}
+	/* do not balance on the way back up */
 
 #ifdef KD_DEBUG
 	/* debug directions */
@@ -1300,52 +1282,57 @@ static struct kdnode *kdtree_insert2(struct kdtree *t, struct kdnode *r,
 #endif
     }
 
-    if (balance && bmode == 0) {
+    if (balance) {
+	int iter, bmode2;
+
+	/* fix any inconsistencies in the (sub-)tree */
+	iter = 0;
+	bmode2 = 0;
 	top = 0;
-	go_back = 0;
 	s[top].n = r;
-	while (s[top].n) {
-	    n = s[top].n;
-	    old_depth = n->depth;
+	while (top >= 0) {
 
-	    /* balance node */
-	    while (kdtree_balance(t, n->child[0], bmode));
-	    while (kdtree_balance(t, n->child[1], bmode));
+	    n = s[top].n;
 
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
-	    while (kdtree_balance(t, n, bmode)) {
-		/*
-		while (kdtree_balance(t, n->child[0], bmode));
-		while (kdtree_balance(t, n->child[1], bmode));
-		*/
-
-		ld = (!n->child[0] ? -1 : n->child[0]->depth);
-		rd = (!n->child[1] ? -1 : n->child[1]->depth);
-		n->depth = MAX(ld, rd) + 1;
+	    /* top-down balancing
+	     * slower but more compact */
+	    if (!bmode2) {
+		while (kdtree_balance(t, n, bmode));
 	    }
-	    if (old_depth != n->depth)
-		go_back = top;
-
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    
-	    dir = (ld < rd);
 
-	    top++;
-	    s[top].n = n->child[dir];
-	}
-
-	top = go_back;
-	while (top) {
-	    top--;
-	    n = s[top].n;
+	    /* go down */
+	    if (n->child[0] && n->child[0]->balance) {
+		dir = 0;
+		top++;
+		s[top].n = n->child[dir];
+	    }
+	    else if (n->child[1] && n->child[1]->balance) {
+		dir = 1;
+		top++;
+		s[top].n = n->child[dir];
+	    }
+	    /* go back up */
+	    else {
 
-	    /* update node depth */
-	    ld = (!n->child[0] ? -1 : n->child[0]->depth);
-	    rd = (!n->child[1] ? -1 : n->child[1]->depth);
-	    n->depth = MAX(ld, rd) + 1;
+		/* bottom-up balancing
+		 * faster but less compact */
+		if (bmode2) {
+		    while (kdtree_balance(t, n, bmode));
+		}
+		top--;
+		if (top >= 0) {
+		    kdtree_update_node(t, s[top].n);
+		}
+		if (!bmode2 && top == 0) {
+		    iter++;
+		    if (iter == 2) {
+			/* the top node has been visited twice, 
+			 * switch from top-down to bottom-up balancing */
+			iter = 0;
+			bmode2 = 1;
+		    }
+		}
+	    }
 	}
     }
 

lib/btree2/kdtree.h

@@ -68,6 +68,7 @@ struct kdnode
 {
     unsigned char dim;          /*!< split dimension of this node */
     unsigned char depth;        /*!< depth at this node */
+    unsigned char balance;      /*!< flag to indicate if balancing is needed */
     double *c;                  /*!< coordinates */
     int uid;                    /*!< unique id of this node */
     struct kdnode *child[2];    /*!< link to children: `[0]` for smaller, `[1]` for larger */