Spatial queries at scale
How BIM software answers “what is here?” for hundreds of thousands of elements: uniform grids, quadtrees and bounding volume hierarchies for window queries, picking rays, clash pairs and nearest neighbours.
The uniform grid
var items = Site.Boxes(20_000);var grid = new GridIndex(5.0);for (int i = 0; i < items.Count; i++) grid.Insert(i, items[i]); var window = new Box(100, 100, 130, 120);int brute = items.Count(b => b.Overlaps(window));var candidates = grid.Query(window);int found = candidates.Count(i => items[i].Overlaps(window)); Console.WriteLine($"brute force: {items.Count} boxes tested, {brute} found");Console.WriteLine($"grid: {grid.CellsTouched(window)} cells opened, {candidates.Count} boxes tested, {found} found");Console.WriteLine($"{grid.CellCount} cells in use, {grid.Entries / (double)items.Count:0.00} cells per element"); sealed class GridIndex{ readonly double size; readonly Dictionary<(int, int), List<int>> cells = new(); public GridIndex(double size) => this.size = size; public int Entries { get; private set; } public int CellCount => cells.Count; (int X0, int Y0, int X1, int Y1) Range(Box b) => ((int)Math.Floor(b.X0 / size), (int)Math.Floor(b.Y0 / size), (int)Math.Floor(b.X1 / size), (int)Math.Floor(b.Y1 / size)); public void Insert(int id, Box b) { var r = Range(b); for (int x = r.X0; x <= r.X1; x++) for (int y = r.Y0; y <= r.Y1; y++) { if (!cells.TryGetValue((x, y), out var list)) cells[(x, y)] = list = new(); list.Add(id); Entries++; } } public int CellsTouched(Box b) { var r = Range(b); return (r.X1 - r.X0 + 1) * (r.Y1 - r.Y0 + 1); } public HashSet<int> Query(Box b) { var found = new HashSet<int>(); // a set removes elements listed in several cells var r = Range(b); for (int x = r.X0; x <= r.X1; x++) for (int y = r.Y0; y <= r.Y1; y++) if (cells.TryGetValue((x, y), out var list)) found.UnionWith(list); return found; }} readonly record struct Box(double X0, double Y0, double X1, double Y1){ public bool Overlaps(Box o) => X0 < o.X1 && o.X0 < X1 && Y0 < o.Y1 && o.Y0 < Y1;} static class Site{ // 20 000 elements between 0.2 m and 2 m across, spread over a 500 x 300 m site; the same every run public static List<Box> Boxes(int n) { uint s = 42; double Next() { s = s * 1664525 + 1013904223; return s / 4294967296.0; } var list = new List<Box>(); for (int i = 0; i < n; i++) { double x = Next() * 500, y = Next() * 300, w = 0.2 + Next() * 1.8, h = 0.2 + Next() * 1.8; list.Add(new Box(x, y, x + w, y + h)); } return list; }}brute force: 20000 boxes tested, 100 found grid: 35 cells opened, 136 boxes tested, 100 found 6046 cells in use, 1.48 cells per element
Trees that adapt: quadtrees and BVHs
var items = Site.Boxes(20_000);var bvh = new Bvh(items);Console.WriteLine($"BVH: {bvh.Nodes} nodes, {bvh.Depth} levels"); var window = new Box(100, 100, 130, 120);var (hits, visited, tested) = bvh.Query(window);Console.WriteLine($"window: visited {visited} nodes, tested {tested} boxes, found {hits}"); // A pick ray in plan: from (0, 150) heading east, slightly northdouble ox = 0, oy = 150, dx = 1, dy = 0.05;var ray = bvh.Ray(ox, oy, dx, dy);int bruteId = Enumerable.Range(0, items.Count).MinBy(i => Bvh.Slab(items[i], ox, oy, dx, dy));Console.WriteLine($"ray: first hit element {ray.Id} at t = {ray.T:0.000}, visited {ray.Visited} nodes, tested {ray.Tested} boxes");Console.WriteLine($"brute force agrees: {bruteId == ray.Id}"); sealed class Bvh{ sealed class Node { public Box B; public Node? L, R; public int[]? Items; } readonly List<Box> items; readonly Node root; public int Nodes { get; private set; } public int Depth { get; private set; } public Bvh(List<Box> items) { this.items = items; root = Build(Enumerable.Range(0, items.Count).ToArray(), 1); } Node Build(int[] ids, int depth) { Nodes++; Depth = Math.Max(Depth, depth); var b = new Box(ids.Min(i => items[i].X0), ids.Min(i => items[i].Y0), ids.Max(i => items[i].X1), ids.Max(i => items[i].Y1)); if (ids.Length <= 4) return new Node { B = b, Items = ids }; bool alongX = b.X1 - b.X0 >= b.Y1 - b.Y0; // split across the longest side var sorted = ids.OrderBy(i => alongX ? items[i].X0 + items[i].X1 : items[i].Y0 + items[i].Y1).ToArray(); int mid = sorted.Length / 2; return new Node { B = b, L = Build(sorted[..mid], depth + 1), R = Build(sorted[mid..], depth + 1) }; } public (int Hits, int Visited, int Tested) Query(Box q) { int hits = 0, visited = 0, tested = 0; var stack = new Stack<Node>(); stack.Push(root); while (stack.Count > 0) { var n = stack.Pop(); visited++; if (!n.B.Overlaps(q)) continue; // the whole branch is outside if (n.Items is null) { stack.Push(n.L!); stack.Push(n.R!); continue; } foreach (int i in n.Items) { tested++; if (items[i].Overlaps(q)) hits++; } } return (hits, visited, tested); } public (int Id, double T, int Visited, int Tested) Ray(double ox, double oy, double dx, double dy) { int best = -1, visited = 0, tested = 0; double bestT = double.PositiveInfinity; var stack = new Stack<Node>(); stack.Push(root); while (stack.Count > 0) { var n = stack.Pop(); visited++; if (Slab(n.B, ox, oy, dx, dy) >= bestT) continue; // cannot beat the best hit so far if (n.Items is null) { // visit the nearer child first so a close hit prunes the other one double tl = Slab(n.L!.B, ox, oy, dx, dy), tr = Slab(n.R!.B, ox, oy, dx, dy); if (tl < tr) { stack.Push(n.R); stack.Push(n.L); } else { stack.Push(n.L); stack.Push(n.R); } continue; } foreach (int i in n.Items) { tested++; double t = Slab(items[i], ox, oy, dx, dy); if (t < bestT) { bestT = t; best = i; } } } return (best, bestT, visited, tested); } // Slab test: where the ray enters the box, or infinity if it misses public static double Slab(Box b, double ox, double oy, double dx, double dy) { double t0 = 0, t1 = double.PositiveInfinity; foreach (var (o, d, lo, hi) in new[] { (ox, dx, b.X0, b.X1), (oy, dy, b.Y0, b.Y1) }) { if (Math.Abs(d) < 1e-15) { if (o < lo || o > hi) return double.PositiveInfinity; continue; } double a = (lo - o) / d, c = (hi - o) / d; t0 = Math.Max(t0, Math.Min(a, c)); t1 = Math.Min(t1, Math.Max(a, c)); } return t0 <= t1 ? t0 : double.PositiveInfinity; }} readonly record struct Box(double X0, double Y0, double X1, double Y1){ public bool Overlaps(Box o) => X0 < o.X1 && o.X0 < X1 && Y0 < o.Y1 && o.Y0 < Y1;} static class Site{ public static List<Box> Boxes(int n) { uint s = 42; double Next() { s = s * 1664525 + 1013904223; return s / 4294967296.0; } var list = new List<Box>(); for (int i = 0; i < n; i++) { double x = Next() * 500, y = Next() * 300, w = 0.2 + Next() * 1.8, h = 0.2 + Next() * 1.8; list.Add(new Box(x, y, x + w, y + h)); } return list; }}BVH: 15423 nodes, 14 levels window: visited 147 nodes, tested 114 boxes, found 100 ray: first hit element 3668 at t = 2.958, visited 47 nodes, tested 7 boxes brute force agrees: True
Clash pairs: the broad phase at scale
var items = Site.Boxes(20_000); long bruteTests = 0, bruteFound = 0;for (int i = 0; i < items.Count; i++) for (int j = i + 1; j < items.Count; j++) { bruteTests++; if (items[i].Overlaps(items[j])) bruteFound++; }Console.WriteLine($"brute force: {bruteTests:N0} pair tests, {bruteFound} overlapping pairs"); const double Size = 2.5;var cells = new Dictionary<(int, int), List<int>>();for (int i = 0; i < items.Count; i++){ var b = items[i]; for (int x = (int)Math.Floor(b.X0 / Size); x <= (int)Math.Floor(b.X1 / Size); x++) for (int y = (int)Math.Floor(b.Y0 / Size); y <= (int)Math.Floor(b.Y1 / Size); y++) { if (!cells.TryGetValue((x, y), out var list)) cells[(x, y)] = list = new(); list.Add(i); }} long gridTests = 0, gridFound = 0;foreach (var ((cx, cy), list) in cells) for (int a = 0; a < list.Count; a++) for (int b = a + 1; b < list.Count; b++) { gridTests++; Box p = items[list[a]], q = items[list[b]]; if (!p.Overlaps(q)) continue; // report the pair only in the cell holding the lower-left corner of the overlap double ox = Math.Max(p.X0, q.X0), oy = Math.Max(p.Y0, q.Y0); if ((int)Math.Floor(ox / Size) == cx && (int)Math.Floor(oy / Size) == cy) gridFound++; }Console.WriteLine($"grid: {gridTests:N0} pair tests, {gridFound} overlapping pairs, {(double)bruteTests / gridTests:0} times fewer tests"); readonly record struct Box(double X0, double Y0, double X1, double Y1){ public bool Overlaps(Box o) => X0 < o.X1 && o.X0 < X1 && Y0 < o.Y1 && o.Y0 < Y1;} static class Site{ public static List<Box> Boxes(int n) { uint s = 42; double Next() { s = s * 1664525 + 1013904223; return s / 4294967296.0; } var list = new List<Box>(); for (int i = 0; i < n; i++) { double x = Next() * 500, y = Next() * 300, w = 0.2 + Next() * 1.8, h = 0.2 + Next() * 1.8; list.Add(new Box(x, y, x + w, y + h)); } return list; }}brute force: 199,990,000 pair tests, 6498 overlapping pairs grid: 35,297 pair tests, 6498 overlapping pairs, 5666 times fewer tests
| Structure | Build | Strengths | Weaknesses |
|---|---|---|---|
Show a hint
Show the solution
var items = Site.Boxes(20_000);var centres = items.Select(b => ((b.X0 + b.X1) / 2, (b.Y0 + b.Y1) / 2)).ToList();const double Size = 5;var grid = new Dictionary<(int, int), List<int>>();for (int i = 0; i < centres.Count; i++){ var key = ((int)Math.Floor(centres[i].Item1 / Size), (int)Math.Floor(centres[i].Item2 / Size)); if (!grid.TryGetValue(key, out var list)) grid[key] = list = new(); list.Add(i);} double qx = 250.3, qy = 150.7;int cx = (int)Math.Floor(qx / Size), cy = (int)Math.Floor(qy / Size);int best = -1, examined = 0;double bestD = double.PositiveInfinity;for (int k = 0; ; k++){ for (int i = -k; i <= k; i++) for (int j = -k; j <= k; j++) { if (Math.Max(Math.Abs(i), Math.Abs(j)) != k) continue; // only the ring itself if (!grid.TryGetValue((cx + i, cy + j), out var list)) continue; foreach (int id in list) { examined++; double d = Math.Sqrt(Math.Pow(centres[id].Item1 - qx, 2) + Math.Pow(centres[id].Item2 - qy, 2)); if (d < bestD) { bestD = d; best = id; } } } if (bestD <= k * Size) break;}int brute = Enumerable.Range(0, centres.Count).MinBy(i => Math.Pow(centres[i].Item1 - qx, 2) + Math.Pow(centres[i].Item2 - qy, 2));Console.WriteLine($"nearest element {best}, {bestD:0.000} m away");Console.WriteLine($"examined {examined} of {centres.Count} centres; brute force agrees: {brute == best}"); readonly record struct Box(double X0, double Y0, double X1, double Y1); static class Site{ public static List<Box> Boxes(int n) { uint s = 42; double Next() { s = s * 1664525 + 1013904223; return s / 4294967296.0; } var list = new List<Box>(); for (int i = 0; i < n; i++) { double x = Next() * 500, y = Next() * 300, w = 0.2 + Next() * 1.8, h = 0.2 + Next() * 1.8; list.Add(new Box(x, y, x + w, y + h)); } return list; }}nearest element 17464, 1.040 m away examined 34 of 20000 centres; brute force agrees: True
Remember
A spatial index stores bounding boxes so a query only tests elements near it, instead of all n. A uniform grid is simple and fast for evenly spread elements; pick a cell a little bigger than a typical element. Quadtrees and octrees split space where elements crowd; BVHs group elements and store each one exactly once. Ray queries on a BVH visit the nearer child first and skip nodes further away than the best hit so far. Grid-based clash search needs a de-duplication rule, such as reporting a pair only in the cell holding the corner of its overlap.