31:T8fd,md`# Omohundro balltree construction algorithms [Balltrees](https://en.wikipedia.org/wiki/Ball_tree) are hierarchical spatial data structures that can be used to answer proximity queries. In that sense, they are like quadtrees, except that they are not a spatial decomposition scheme, but rather rely on spatial aggregation. Unlike [quadtrees](https://en.wikipedia.org/wiki/Quadtree), one point in space can be inside many nodes at the same level of a balltree, which makes them less adequate for indexing points. On the other hand, when indexing objects with extent, i.e., non-zero size, they can overlap more than one leaf quatree node, whereas balltrees need only one leaf pointing to each object. In this, they are perhaps more related to [R-trees](https://en.wikipedia.org/wiki/R-tree). In this notebook, we are interested in indexing collections of balls of various sizes in 2D. When balls have very small radii, they can be regarded as points. We have implemented four of the [Five balltree construction algorithms](http://www.icsi.berkeley.edu/ftp/global/pub/techreports/1989/tr-89-063.pdf), as described in that seminal paper by Stephen Omohundro: 1. *KD construction*. A method reminiscent of the one used for building offline [k-d trees](https://en.wikipedia.org/wiki/K-d_tree). The collection of objects is recursively split by the median along the axis of maximum dispersion. 2. *Online construction*. A method for building balltrees incrementally, supporting both insertion and deletion. 3. *Bottom up construction*. A method whereby the tree is built by recursively pairing the balls which would yield the smallest enclosing ball. The brute-force algorithm is quite expensive - ${tex`O(n^3)`} - and thus Omohundro also proposes a ... 4. *Clever bottom up construction* algorithm, where a priority queue and an auxiliary Online balltree are used to reduce the complexity of the brute-force construction algorithm. Note that there are two other algorithms proposed by Omohundro: The *top-down* and the *cheap online* aproaches. These, however, tend to yield less efficient trees. Incidentally, the main measure used to gauge the efficiency of ball trees is that of total area (or volume for dimensions higher than 2), i.e., sum of the areas taken by every node in the tree. `32:T177b,// // Ball base class // class Ball { constructor(x, y, r) { [this.x, this.y, this.r] = [x, y, r]; } // Creates a copy of other ball static from(other) { return new Ball(other.x, other.y, other.r); } // Returns true iff this and other are "equal", i.e., same center and radius equals(other) { const dx = this.x - other.x; const dy = this.y - other.y; const dr = this.r - other.r; return dx * dx == 0 && dy * dy == 0 && dr * dr == 0; } // Returns true iff this other ⊆ this. encloses(other) { const dr = this.r - other.r + 1e-6; // const dx = this.x - other.x; const dy = this.y - other.y; return dr >= 0 && dr * dr > dx * dx + dy * dy; } // returns true iff this ball contains point p containsPoint(p) { const dx = this.x - p.x; const dy = this.y - p.y; return dx * dx + dy * dy <= this.r * this.r; } // Returns true iff this ball intersects other intersects(other) { const dx = this.x - other.x; const dy = this.y - other.y; return Math.sqrt(dx * dx + dy * dy) < this.r + other.r; } // Returns true iff this b ⊆ this for all b ∈ others. enclosesAll(others) { for (let b of others) if (!this.encloses(b)) return false; return true; } // returns true iff the smallest enclosing ball of a and b is neither a nor b static isBasis2(a, b) { return !a.encloses(b) && !b.encloses(a); } // returns smallest enclosing ball containing this and other union(other) { if (other.encloses(this)) return other; if (this.encloses(other)) return this; const { x: x1, y: y1, r: r1 } = this; const { x: x2, y: y2, r: r2 } = other; const x21 = x2 - x1, y21 = y2 - y1, r21 = r2 - r1; const l = Math.sqrt(x21 * x21 + y21 * y21); return new Ball( (x1 + x2 + (x21 / l) * r21) / 2, (y1 + y2 + (y21 / l) * r21) / 2, (l + r1 + r2) / 2 ); } // Returns this ball expanded so that it includes at least one point of other expand(other) { const { x: x1, y: y1, r: r1 } = this; const { x: x2, y: y2, r: r2 } = other; const x21 = x2 - x1, y21 = y2 - y1, r21 = r2 - r1; const l = Math.sqrt(x21 * x21 + y21 * y21); if (l <= this.r + other.r) return this; const R = (l - r2 + r1) / 2; return new Ball(x1 + (x21 * (R - r1)) / l, y1 + (y21 * (R - r1)) / l, R); } // Returns true iff the set {a, b, c} forms a basis. static isBasis3(a, b, c) { return ( Ball.isBasis2(a, b) && !a.union(b).encloses(c) && Ball.isBasis2(a, c) && !a.union(c).encloses(b) && Ball.isBasis2(b, c) && !b.union(c).encloses(a) ); } // If three balls form a 'basis' (see isBasis3), computes the // smallest enclosing circle of these static encloseBasis3(a, b, c) { const { x: x1, y: y1, r: r1 } = a; const { x: x2, y: y2, r: r2 } = b; const { x: x3, y: y3, r: r3 } = c; const a2 = x1 - x2; const a3 = x1 - x3; const b2 = y1 - y2; const b3 = y1 - y3; const c2 = r2 - r1; const c3 = r3 - r1; const d1 = x1 * x1 + y1 * y1 - r1 * r1; const d2 = d1 - x2 * x2 - y2 * y2 + r2 * r2; const d3 = d1 - x3 * x3 - y3 * y3 + r3 * r3; const ab = a3 * b2 - a2 * b3; const xa = (b2 * d3 - b3 * d2) / (ab * 2) - x1; const xb = (b3 * c2 - b2 * c3) / ab; const ya = (a3 * d2 - a2 * d3) / (ab * 2) - y1; const yb = (a2 * c3 - a3 * c2) / ab; const A = xb * xb + yb * yb - 1; const B = 2 * (r1 + xa * xb + ya * yb); const C = xa * xa + ya * ya - r1 * r1; /** ** This deviates from the implementation in d3.packEnclose(), which ** throws an error for ** ** d3.packEnclose ([ ** { x: 14.5, y: 48.5, r: 7.585 }, ** { x: 9.5, y: 79.5, r: 2.585 }, ** { x: 15.5, y: 73.5, r: 8.585 } ** ]) ** ** (see https://github.com/d3/d3-hierarchy/issues/188) **/ const r = -(Math.abs(A) > 1e-10 ? (B + Math.sqrt(B * B - 4 * A * C)) / (2 * A) : C / B); return new Ball(x1 + xa + xb * r, y1 + ya + yb * r, r); } // Given a basis B and a circle p ⊈ B, returns the new basis Bʹ. static extendBasis(B, p) { var i, j; if (p.enclosesAll(B)) return [p]; // If we get here then B must have at least one element. for (i = 0; i < B.length; ++i) { if (!p.encloses(B[i]) && p.union(B[i]).enclosesAll(B)) { return [B[i], p]; } } // If we get here then B must have at least two elements. for (i = 0; i < B.length - 1; ++i) { for (j = i + 1; j < B.length; ++j) { if ( !B[i].union(B[j]).encloses(p) && !B[i].union(p).encloses(B[j]) && !B[j].union(p).encloses(B[i]) && Ball.encloseBasis3(B[i], B[j], p).enclosesAll(B) ) { return [B[i], B[j], p]; } } } // If we get here then something is very wrong. throw "Something Weird happened"; } // Returns a ball enclosing all balls in L static enclose(L) { // let { x, y, r } = d3.packEnclose(L); // return new Ball(x, y, r); const encloseBasis = (B) => { switch (B.length) { case 1: return B[0]; case 2: return B[0].union(B[1]); case 3: return Ball.encloseBasis3(B[0], B[1], B[2]); } }; var i = 0, n = d3.shuffle((L = L.slice())).length, B = [], p, e; while (i < n) { p = L[i]; if (e && e.encloses(p)) ++i; else (e = encloseBasis((B = Ball.extendBasis(B, p)))), (i = 0); } return e; } // Returns the distance from a point to this ball dist(p) { return Math.max( 0, Math.sqrt((this.x - p.x) ** 2 + (this.y - p.y) ** 2) - this.r ); } // Area of this ball get area() { return this.r * this.r * Math.PI; } draw(ctx, fill = true, stroke = true) { ctx.beginPath(); ctx.arc(this.x, this.y, this.r, 0, Math.PI * 2); if (stroke) ctx.stroke(); if (fill) ctx.fill(); } }33:T1eda,// // Balltree base class // class BT { // Builds a balltree contaning a BTnode (or an empty balltree) constructor(node = null) { this.node = node; } // Height of the tree height() { if (!this.node) return -1; if (this.node.leaf()) return 0; return Math.max(this.node.left.height() + 1, this.node.right.height() + 1); } // Visits the tree in inorder, yielding (generating) all nodes *inorder() { if (this.node) { if (this.node.left) yield* this.node.left.inorder(); yield this.node; if (this.node.right) yield* this.node.right.inorder(); } } // Visits the tree generating all nodes at level l *level(l) { if (l == 0) yield this.node; else { if (this.node.left) yield* this.node.left.level(l - 1); if (this.node.right) yield* this.node.right.level(l - 1); } } // Similar to level, but includes leaf nodes if at level shallower than l *approximation(l) { if (l == 0 || this.node.leaf()) yield this.node; else { if (this.node.left) yield* this.node.left.approximation(l - 1); if (this.node.right) yield* this.node.right.approximation(l - 1); } } // Generates all leaf nodes of this tree *leafBalls() { if (this.node.leaf()) yield this.node.entry; else { yield* this.node.left.leafBalls(); yield* this.node.right.leafBalls(); } } // Generates all leaf notes within distance d of point p *within(p, d) { let self = this; self.withinSteps = 0; function* traverse(node) { let dnode = node.entry.dist(p); self.withinSteps++; if (dnode <= d) { if (node.leaf()) yield node; else { let dleft = node.left.node.entry.dist(p); let dright = node.right.node.entry.dist(p); if (dleft < dright) { yield* traverse(node.left.node); yield* traverse(node.right.node); } else { yield* traverse(node.left.node); yield* traverse(node.right.node); } } } } if (this.node) yield* traverse(this.node); } // Reduces the bounding balls from internal nodes to those // enclosing only the leaves computeEnclosingCircles() { if (!this.node.leaf()) { this.node.entry = Ball.enclose([...this.leafBalls()]); //Ball.from(d3.packEnclose([...this.leafBalls()])); this.node.left.computeEnclosingCircles(); this.node.right.computeEnclosingCircles(); } return this; } // Total area of the tree area() { if (this.node.leaf()) return this.node.entry.area; return ( this.node.entry.area + this.node.left.area() + this.node.right.area() ); } // Draws all balls using context ctx draw(ctx) { for (let node of this.inorder()) { node.entry.draw(ctx); } } // Finds a leaf matching position and radius with the given ball // returns a pair {node, ancestors} where node is // the leaf node with the found ball, and ancestor is an array of // tree objects [T0, T1, ... Tn], such that Ti+1 is the parent of Ti, // Tn is the tree containing node find(ball) { this.findSteps = 0; this.visited = []; let traverse = (ancestors) => { let n = ancestors.length; let node = ancestors[n - 1].node; this.visited.push(node); this.findSteps++; if (node.leaf()) { if (node.entry.equals(ball)) return { node, ancestors }; return null; } let result = null; if (node.entry.containsPoint(ball)) { let dleft = node.left.node.entry.dist(ball); let dright = node.right.node.entry.dist(ball); if (dleft < dright) { ancestors.push(node.left); result = traverse(ancestors); if (result) return result; ancestors.pop(); ancestors.push(node.right); result = traverse(ancestors); } else { ancestors.push(node.right); result = traverse(ancestors); if (result) return result; ancestors.pop(); ancestors.push(node.left); result = traverse(ancestors); } if (result) return result; ancestors.pop(); } return null; }; if (this.node) return traverse([this]); } // Returns a leaf ball containing p, or null, if none exists intersecting(p) { this.intersectingSteps = 0; this.visited = []; let traverse = (node) => { this.intersectingSteps++; this.visited.push(node); if (node.leaf()) { if (node.entry.containsPoint(p)) return node.entry; return null; } if (!node.entry.containsPoint(p)) return null; let dleft = node.left.node.entry.dist(p); let dright = node.right.node.entry.dist(p); if (dleft < dright) { let result = traverse(node.left.node); if (result) return result; return traverse(node.right.node); } let result = traverse(node.right.node); if (result) return result; return traverse(node.left.node); }; return traverse(this.node); } // Returns a list of balls enclosed in the given ball findInCircle(ball) { this.findInCircleSteps = 0; this.visited = []; let result = []; let traverse = (node) => { this.findInCircleSteps++; this.visited.push(node); if (node.leaf()) { if (ball.encloses(node.entry)) result.push(node.entry); } else { if (node.left.node.entry.intersects(ball)) traverse(node.left.node); if (node.right.node.entry.intersects(ball)) traverse(node.right.node); } }; if (this.node.entry.intersects(ball)) traverse(this.node); return result; } // Returns the ball in a leaf node whose center is closest to point p. // Uses branch and bound to speed up the search. closest(p) { let best = null; let bestDist = Number.POSITIVE_INFINITY; this.closestSteps = 0; this.visited = []; let traverse = (node) => { if (bestDist == 0) return; this.closestSteps++; this.visited.push(node); if (node.leaf()) { let d = node.entry.dist(p); if (d < bestDist) [best, bestDist] = [node.entry, d]; return; } let dleft = node.left.node.entry.dist(p); let dright = node.right.node.entry.dist(p); if (dleft < dright) { if (dleft < bestDist) traverse(node.left.node); if (dright < bestDist) traverse(node.right.node); } else { if (dright < bestDist) traverse(node.right.node); if (dleft < bestDist) traverse(node.left.node); } }; if (this.node) traverse(this.node); return best; } // Returns a leaf node which would result in the smallest union with the given ball. // Uses branch and bound to speed the search smallestUnionLeaf(ball) { let best = null; let bestUnionArea = Number.POSITIVE_INFINITY; this.smallestUnionSteps = 0; this.visited = []; let traverse = (node) => { this.smallestUnionSteps++; this.visited.push(node); let expansion = ball.expand(node.entry).area; if (expansion > bestUnionArea) return; if (node.leaf()) { let union = node.entry.union(ball); let area = union.area; if (area < bestUnionArea) { best = node; bestUnionArea = area; } return; } let eleft = ball.expand(node.left.node.entry); let eright = ball.expand(node.right.node.entry); if (eleft.area < eright.area) { if (eleft.area < bestUnionArea) traverse(node.left.node); if (eright.area < bestUnionArea) traverse(node.right.node); } else { if (eright.area < bestUnionArea) traverse(node.right.node); if (eleft.area < bestUnionArea) traverse(node.left.node); } }; if (this.node) traverse(this.node); return best; } }34:Tef6,class OnlineBT extends BT { constructor(...args) { super(...args); } // Finds a leaf having the given ball as entry (must be the same object!). // Returns a pair {node, ancestors} where node is // the leaf node with the found ball, and ancestor is an array of // tree objects [T0, T1, ... Tn], such that Ti+1 is the parent of Ti, // Tn is the tree containing node find(ball) { this.findSteps = 0; this.visited = []; let traverse = ancestors => { let n = ancestors.length; let node = ancestors[n - 1].node; this.visited.push(node); this.findSteps++; if (node.leaf()) { if (node.entry == ball) return { node, ancestors }; return null; } let result = null; if (node.entry.containsPoint(ball)) { let dleft = node.left.node.entry.dist(ball); let dright = node.right.node.entry.dist(ball); if (dleft < dright) { ancestors.push(node.left); result = traverse(ancestors); if (result) return result; ancestors.pop(); ancestors.push(node.right); result = traverse(ancestors); } else { ancestors.push(node.right); result = traverse(ancestors); if (result) return result; ancestors.pop(); ancestors.push(node.left); result = traverse(ancestors); } if (result) return result; ancestors.pop(); } return null; }; if (this.node) return traverse([this]); } remove(ball) { let findResult = this.find(ball); if (!findResult) { throw "Trying to remove inexistant ball"; } let { node, ancestors } = findResult; if (ancestors.length > 1) { ancestors.pop(); // Pop tree containing ball let parent = ancestors[ancestors.length - 1]; if (node == parent.node.left.node) { parent.node = parent.node.right.node; } else { parent.node = parent.node.left.node; } ancestors.pop(); // Parent tree now has a leaf node // Update ancestors of parent while (ancestors.length > 0) { let tree = ancestors.pop(); tree.node.entry = tree.node.left.node.entry.union( tree.node.right.node.entry ); } } else { // Root node was removed this.node = null; } } insert(ball) { if (!this.node) { this.node = new BTNode(ball); return; } let PQNode = (tree, ancestorExpansion, parentpq) => { let union = tree.node.entry.union(ball); let cost = union.area + ancestorExpansion; let aexp = ancestorExpansion; return { tree, union, cost, aexp, parentpq }; }; let pq = new BinaryHeap(pqnode => pqnode.aexp); let best = PQNode(this, 0, null); pq.push(best); while (pq.size() > 0) { let pqnode = pq.pop(); if (pqnode.aexp >= best.cost) break; let e = pqnode.aexp + pqnode.cost - pqnode.tree.node.entry.area; if (pqnode.tree.node.left) { let left = PQNode(pqnode.tree.node.left, e, pqnode); if (left.cost < best.cost) best = left; pq.push(left, e); } if (pqnode.tree.node.right) { let right = PQNode(pqnode.tree.node.right, e, pqnode); if (right.cost < best.cost) best = right; pq.push(right, e); } } // Make ball a sibling of best node best.tree.node = new BTNode( best.tree.node.entry.union(ball), //best.union, new OnlineBT(best.tree.node), new OnlineBT(new BTNode(ball)) ); // Update ancestor bounding balls let pqnode = best; while (pqnode.parentpq) { pqnode = pqnode.parentpq; let node = pqnode.tree.node; let union = node.left.node.entry.union(node.right.node.entry); node.entry = union; } } }35:T7ca,function bottomUp(balls) { let trees = balls.map(b => new BT(new BTNode(b))); let PQNode = (i, j, unionArea) => ({ i, j, unionArea }); let pq = new BinaryHeap(pqnode => pqnode.unionArea); let auxtree = new OnlineBT(); // Build auxliary tree and priority queue for (let i = 0; i < trees.length; i++) { let best = null; let bestArea = Number.POSITIVE_INFINITY; for (let j = 0; j < trees.length; j++) { if (i == j) continue; let unionArea = trees[i].node.entry.union(trees[j].node.entry).area; if (unionArea < bestArea) { bestArea = unionArea; best = PQNode(i, j, bestArea); } } pq.push(best); trees[i].node.entry.index = i; auxtree.insert(trees[i].node.entry); } // insert pairs from the priority queue pq let iteration = 0; while (pq.size() > 1) { iteration++; let pqnode = pq.pop(); let pqtop = pq.content[0]; let { i, unionArea } = pqnode; // Check whether this ball has already been inserted if (!trees[i]) continue; // Recompute the best mate for tree i let ball = trees[i].node.entry; auxtree.remove(ball); let mate = auxtree.smallestUnionLeaf(ball).entry; let union = Ball.from(ball.union(mate)); // Make a copy if (union.area <= pqtop.unionArea) { let parent = new BT( new BTNode(union, trees[ball.index], trees[mate.index]) ); trees[ball.index] = null; trees[mate.index] = null; auxtree.remove(mate); union.index = trees.length; trees.push(parent); let bestParentMate = auxtree.smallestUnionLeaf(union); if (!bestParentMate) break; auxtree.insert(union); pq.push( PQNode( union.index, bestParentMate.entry.index, union.union(bestParentMate.entry).area ) ); } else { auxtree.insert(ball); pq.push(PQNode(i, mate.index, union.area)); } } let pqnode = pq.pop(); return trees[trees.length - 1]; }36:T41e,function kdconstruct(balls) { let trees = balls.map(b => new BT(new BTNode(b))); const bbox = trees => { let min = { x: Number.POSITIVE_INFINITY, y: Number.POSITIVE_INFINITY }; let max = { x: Number.NEGATIVE_INFINITY, y: Number.NEGATIVE_INFINITY }; for (let tree of trees) { let ball = tree.node.entry; min.x = Math.min(ball.x, min.x); min.y = Math.min(ball.y, min.y); max.x = Math.max(ball.x, max.x); max.y = Math.max(ball.y, max.y); } return { min, max }; }; const build = trees => { if (trees.length == 1) return trees[0]; let { min, max } = bbox(trees); if (max.x - min.x > max.y - min.y) { medianSort(trees, tree => tree.node.entry.x); } else { medianSort(trees, tree => tree.node.entry.y); } let k = trees.length >> 1; let left = build(trees.slice(0, k)); let right = build(trees.slice(k)); return new BT( new BTNode(left.node.entry.union(right.node.entry), left, right) ); }; let result = build(trees); return result; }37:T8b3,function quadratic(balls) { let trees = balls.map((b) => new BT(new BTNode(b))); const treeSort = (trees) => { // random order // return d3.shuffle(trees); // sort in order of size trees.sort((a, b) => a.node.entry.r - b.node.entry.r); return trees; // sort in order of position w.r.t. largest dispersion axis let min = { x: Number.POSITIVE_INFINITY, y: Number.POSITIVE_INFINITY }; let max = { x: Number.NEGATIVE_INFINITY, y: Number.NEGATIVE_INFINITY }; for (let t of trees) { let { x, y } = t.node.entry; if (x < min.x) min.x = x; if (y < min.y) min.y = y; if (x > max.x) max.x = x; if (y > max.y) max.y = y; } if (max.x - min.x > max.y - min.y) trees.sort( (a, b) => a.node.entry.x - a.node.entry.r - b.node.entry.x + b.node.entry.r ); else trees.sort( (a, b) => a.node.entry.y - a.node.entry.r - b.node.entry.y + b.node.entry.r ); }; while (trees.length > 1) { treeSort(trees); let before = trees.length; let newTrees = [], remainder = []; for (let i = 0; i < trees.length; i++) { if (!trees[i]) continue; let best = null; let bestArea = Number.POSITIVE_INFINITY; for (let j = i + 1; j < trees.length; j++) { if (!trees[j]) continue; let union = trees[i].node.entry.union(trees[j].node.entry); if (union.area < bestArea) { best = [j, union]; bestArea = union.area; } } if (best) { let [j, union] = best; newTrees.push(new BT(new BTNode(union, trees[i], trees[j]))); trees[j] = trees[i] = null; } else { remainder.push(trees[i]); } } // Use only about half of the new trees let n = Math.min(newTrees.length, 1 + Math.trunc(newTrees.length * 0.5)); newTrees.sort((a, b) => a.node.entry.r - b.node.entry.r); // sort in order of size trees = newTrees.slice(0, n); for (let i = n; i < newTrees.length; i++) { let t = newTrees[i]; trees.push(t.node.left, t.node.right); } for (let t of remainder) trees.push(t); if (trees.length >= before) throw "too slow"; } return trees[0]; }38:Tdc9,function showTree(tree, name = "") { let title = md`### ${name}`; let level = html``; let label = html` Highlight level: `; let area = html` Total area: ${tree.area()}`; let cost = html` Average find cost: ${averageFindCost(tree)}`; let closestSteps = html``; let closestInfo = html` Steps to locate closest leaf to query point: `; let intersectingSteps = html``; let intersectingInfo = html` Steps to locate leaf containing query point: `; let smallestUnionSteps = html``; let smallestUnionInfo = html` Steps to locate leaf that yields smallest union with point: `; closestInfo.append(closestSteps); intersectingInfo.append(intersectingSteps); smallestUnionInfo.append(smallestUnionSteps); let width = 600; let canvas = DOM.canvas(width, width); let ctx = canvas.getContext("2d"); let mat2d = glmatrix.mat2d; let transform = mat2d.fromScaling([], [width, width]); mat2d.translate(transform, transform, [0.25, 0.25]); mat2d.scale(transform, transform, [0.5, 0.5]); let inverse = mat2d.invert([], transform); let queryBall = new Ball(queryPoint.x, queryPoint.y, 0.0); const queryRadius = 0.3; let refresh = () => { ctx.save(); ctx.clearRect(0, 0, width, width); ctx.lineWidth = 1 / width; ctx.fillStyle = "#0000000a"; ctx.setTransform(...transform); tree.draw(ctx); let l = level.value; ctx.strokeStyle = "red"; ctx.fillStyle = "#ff000020"; ctx.lineWidth = 3 / width; for (let node of tree.level(l)) node.entry.draw(ctx); let closest = tree.closest(queryPoint); closestSteps.innerHTML = tree.closestSteps; let intersecting = tree.intersecting(queryPoint); intersectingSteps.innerHTML = (intersecting ? "(successful) " : "(unsuccessful) ") + tree.intersectingSteps; let smallestUnionNode = tree.smallestUnionLeaf(queryBall); let smallestUnion = smallestUnionNode.entry.union(queryBall); smallestUnionSteps.innerHTML = tree.smallestUnionSteps; ctx.fillStyle = "#0000ff70"; ctx.strokeStyle = "blue"; closest.draw(ctx); ctx.fillStyle = "#ffffff70"; if (intersecting) intersecting.draw(ctx); ctx.setLineDash([3 / width, 3 / width]); ctx.strokeStyle = "green"; smallestUnionNode.entry.draw(ctx); smallestUnion.draw(ctx, false, true); ctx.setLineDash([]); ctx.strokeStyle = "blue"; ctx.beginPath(); ctx.moveTo(queryPoint.x - 0.02, queryPoint.y); ctx.lineTo(queryPoint.x + 0.02, queryPoint.y); ctx.moveTo(queryPoint.x, queryPoint.y - 0.02); ctx.lineTo(queryPoint.x, queryPoint.y + 0.02); ctx.stroke(); ctx.setLineDash([5 / width, 5 / width]); ctx.fillStyle = "#00000000"; new Ball(queryBall.x, queryBall.y, queryRadius).draw(ctx); ctx.strokeStyle = "yellow"; ctx.setLineDash([]); for (let node of tree.within(queryPoint, queryRadius)) { node.entry.draw(ctx); } ctx.restore(); }; canvas.onclick = (e) => { let [x, y] = vec2.transformMat2d([], [e.offsetX, e.offsetY], inverse); mutable queryPoint = { x, y }; refresh(); }; refresh(); level.oninput = () => refresh(); return html`

${title}${label}${level}
${area}
${cost}
${closestInfo}
${intersectingInfo}
${smallestUnionInfo}
${canvas}

`; }39:[] 2f:{"title":"Omohundro balltree construction algorithms","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":418,"value":"md`## Examples\nBelow we show example ball trees built with the implemented algorithms. You can configure the size of the input collection of balls and whether to build standard balltrees, where the cirle of each internal node bounds the only the balls or the two immediate descendants, or whether to make bounding circles enclose all descendant leaves. Clicking on each diagram changes a query point ${tex`q`}, shown as a blue cross. Three queries are supported: \n 1. Find closest leaf node to ${tex`q`} - the leaf node is shown as a blue circle;\n 2. Find leaf node intersecting ${tex`q`} - the leaf node, if found, is shown filled with a lighter blue;\n 3. Find leaf node ${tex`N`} which would yield the smallest ball containing both ${tex`N`}and ${tex`q`} - ${tex`N`} is shown with a lighter fill and the resulting enclosing ball is shown with a green dashed outline.\n\nFor all queries, the total number of steps (node visits) required to reach the answer is displayed.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":550,"value":"viewof config = combo({\n nballs: number({\n title: \"Number of balls\",\n min: 2,\n max: 1000,\n step: 1,\n value: 100\n }),\n radiusRange: rangeSlider({\n title: \"Radius range\",\n min: 0.01,\n max: 0.3,\n value: [0.01, 0.03],\n step: 0.01\n }),\n encloseLeaves: checkbox({\n title: \"Bounding circle enclosing leaves\",\n options: [\"y\"],\n value: \"y\"\n })\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":485,"value":"showTree(testKDConstruct, \"KD construction algorithm\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":411,"value":"showTree(testOnlineBT, \"Online algorithm\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":753,"value":"viewof showBruteForce = checkbox({\n title: \"Show brute force Bottom Up tree\",\n description:\n \"Normally, leave this unchecked, since it should be identical to the clever algorithm.\",\n options: [\"y\"]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":253,"value":"showBruteForce\n ? showTree(testBottomUpBF, \"Bottom up (brute force) algorithm\")\n : md`### Bottom up (brute force) algorithm omitted`","pinned":false,"hidden":false,"mode":"ojs"},{"id":732,"value":"showTree(testBottomUp, \"Bottom up (clever) algorithm\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":870,"value":"showTree(testQuadratic, \"Quadratic algorithm\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":10,"value":"md`## Main data structures\n\nEnclosing ball algorithm borrowed from [d3.packenclose](https://observablehq.com/@d3/d3-packenclose).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":131,"value":"//\n// Balltree node base class\n//\nclass BTNode {\n \n // Entry is a Ball, while left and right are BT's\n constructor (entry,left=null,right=null) {\n Object.assign(this,{entry,left,right});\n } \n \n // Returns true iff this is a leaf node\n leaf () {\n return !(this.left || this.right)\n }\n \n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":133,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":586,"value":"// Locates every leaf of the tree starting from the root and returns\n// the average number of steps required\nfunction averageFindCost(tree) {\n let totalSteps = 0;\n let n = 0;\n for (let ball of tree.leafBalls()) {\n n++;\n let result = tree.find(ball);\n console.assert(result);\n totalSteps += tree.findSteps;\n if (!result) throw \"Something Wrong with find\";\n }\n return totalSteps / n;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":84,"value":"md`## Random Ball generator`","pinned":false,"hidden":false,"mode":"ojs"},{"id":86,"value":"function* balls(n, minRadius = 0.01, maxRadius = 0.01) {\n while (n--)\n yield new Ball(\n Math.random(),\n Math.random(),\n Math.random() * (maxRadius - minRadius) + minRadius\n );\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":4,"value":"md`## Online algorithm`","pinned":false,"hidden":false,"mode":"ojs"},{"id":88,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":343,"value":"testOnlineBT = {\n let t0 = performance.now();\n let tree = new OnlineBT();\n for (let ball of testData) {\n tree.insert(ball);\n }\n if (config.encloseLeaves) tree.computeEnclosingCircles();\n let t1 = performance.now();\n mutable times[\"online\"] = t1 - t0;\n return tree;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":378,"value":"md`## Bottom-up algorithm (brute force)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":380,"value":"function bottomUpBF(balls) {\n let trees = balls.map(b => new BT(new BTNode(b)));\n while (trees.length > 1) {\n let best = null,\n bestArea = Number.POSITIVE_INFINITY;\n for (let i = 0; i < trees.length; i++) {\n for (let j = i + 1; j < trees.length; j++) {\n let union = trees[i].node.entry.union(trees[j].node.entry);\n if (union.area < bestArea) {\n best = [i, j, union];\n bestArea = union.area;\n }\n }\n }\n let [i, j, union] = best;\n trees[i] = new BT(new BTNode(union, trees[i], trees[j]));\n trees.splice(j, 1);\n //return { trees, best };\n }\n return trees[0];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":736,"value":"testBottomUpBF = {\n if (showBruteForce) {\n let tree = bottomUpBF(testData);\n if (config.encloseLeaves) tree.computeEnclosingCircles();\n return tree;\n }\n return null;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":659,"value":"md`## Bottom-up algorithm (clever)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":661,"value":"$35","pinned":true,"hidden":false,"mode":"ojs"},{"id":383,"value":"testBottomUp = {\n let t0 = performance.now();\n let tree = bottomUp(testData);\n if (config.encloseLeaves) tree.computeEnclosingCircles();\n let t1 = performance.now();\n mutable times[\"bottom up\"] = t1 - t0;\n return tree;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":429,"value":"md`## KD construction algorithm`","pinned":false,"hidden":false,"mode":"ojs"},{"id":433,"value":"$36","pinned":true,"hidden":false,"mode":"ojs"},{"id":481,"value":"testKDConstruct = {\n let t0 = performance.now();\n let result = kdconstruct(testData);\n if (config.encloseLeaves) result.computeEnclosingCircles();\n let t1 = performance.now();\n mutable times[\"KD\"] = t1 - t0;\n return result;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":434,"value":"// Moves the median of A to its proper place in the array, that is,\n// to position A.length/2. Function f can be used to establish the\n// sort order, i.e., considers f(A[i]) rather than A[i].\n// Returns the sorted array\nfunction medianSort(A, f = x => x) {\n // Partitions A [i ... i+n-1] using A[i] as pivot;\n // Returns the number of elements < pivot;\n const partition = (i, n) => {\n let x = f(A[i]);\n let m = 0;\n for (let j = i + 1; j < i + n; j++) {\n if (f(A[j]) < x) {\n m++;\n [A[j], A[i + m]] = [A[i + m], A[j]];\n }\n }\n [A[i], A[i + m]] = [A[i + m], A[i]];\n return m;\n };\n\n // Recursively find the median using the quicksort rationale\n const select = (k, i, n) => {\n let j = ~~(Math.random() * n) + i;\n [A[i], A[j]] = [A[j], A[i]];\n let m = partition(i, n);\n if (m == k) return m;\n if (m < k) return m + 1 + select(k - m - 1, i + m + 1, n - m - 1);\n return select(k, i, m);\n };\n\n select(A.length >> 1, 0, A.length);\n return A;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":857,"value":"md`## Quadratic algorithm`","pinned":false,"hidden":false,"mode":"ojs"},{"id":860,"value":"$37","pinned":true,"hidden":false,"mode":"ojs"},{"id":864,"value":"testQuadratic = {\n let t0 = performance.now();\n let tree = quadratic(testData);\n if (config.encloseLeaves) tree.computeEnclosingCircles();\n let t1 = performance.now();\n mutable times[\"quadratic\"] = t1 - t0;\n return tree;\n}\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":886,"value":"mutable times = ({})","pinned":true,"hidden":false,"mode":"ojs"},{"id":502,"value":"md`## Array of balls used for the examples`","pinned":false,"hidden":false,"mode":"ojs"},{"id":163,"value":"testData = [...balls(config.nballs, ...config.radiusRange)]","pinned":true,"hidden":false,"mode":"ojs"},{"id":497,"value":"md`## Component to display ball trees`","pinned":false,"hidden":false,"mode":"ojs"},{"id":580,"value":"mutable queryPoint = ({ x: 0, y: 0 })","pinned":true,"hidden":false,"mode":"ojs"},{"id":247,"value":"$38","pinned":true,"hidden":false,"mode":"ojs"},{"id":112,"value":"md`## Dependencies`","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"glmatrix = import('https://unpkg.com/gl-matrix@3.3.0/esm/index.js?module')","pinned":true,"hidden":false,"mode":"ojs"},{"id":118,"value":"vec2 = glmatrix.vec2","pinned":true,"hidden":false,"mode":"ojs"},{"id":93,"value":"import { BinaryHeap } from \"@esperanc/2d-geometry-utils\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":416,"value":"import { number, checkbox } from \"@jashkenas/inputs\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":548,"value":"import { rangeSlider } from \"@mootari/range-slider\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":553,"value":"import { combo } from \"@esperanc/aggregated-inputs\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":168,"value":"d3 = require(\"d3@6\")","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W39","stdlib":"1","titling":"auto"}