Spatial Indexes: k-d Trees and Quadtrees
Answer range and nearest-neighbor queries over many points without scanning them all, with k-d trees, point-region quadtrees and packed R-trees, and know when each one fails.
Many points, many questions
A delivery app keeps the positions of 100,000 drivers. Every few seconds it asks two kinds of question. Which drivers are inside this rectangle on the map? That is a range query. Which driver is closest to this customer? That is a nearest-neighbor query. A game asks the same questions about units, a CAD program about vertices under the mouse, a GIS about buildings inside a district.
The direct answer is a loop over every point. That costs n distance computations per query, so q queries cost n * q. With 100,000 points and 10,000 queries a second, that is a billion distance computations a second, almost all of them for points nowhere near the query.
A spatial index is a data structure that groups nearby points together so that a query can skip whole groups at once. This lesson builds three: the k-d tree, which splits the points at medians; the point-region quadtree, which splits space into equal squares; and the R-tree, which groups rectangles under bounding boxes and is what databases use. For each we count how much work a query does, test it against brute force, and look for the input that makes it slow.
Coordinates are integers, as in the rest of the course from the robustness lesson on: struct Point { long long x, y; }, with |x|, |y| <= 10^9. Distances are compared as squared distances, which are exact in long long under that bound, so no square root and no rounding ever decides which point is nearer.
Halving the plane: the idea of a k-d tree
Sort the points by x and take the median. Draw a vertical line through it. Half of the points are on the left of the line, half on the right. Now take each half, sort it by y, take its median and draw a horizontal line through it, but only across that half. Keep alternating between x and y until every piece holds one point.
The result is a binary tree. Each node is a point and a splitting line; its left subtree holds the points on one side of the line, its right subtree the points on the other. Because each split is at a median, the two subtrees have the same size, give or take one, and the depth is about log2(n). The name comes from Jon Bentley's 1975 paper, where "k-d" means k-dimensional; here k is 2.
A worked example with seven points: (1,6) (3,2) (4,9) (6,5) (8,1) (9,7) (11,4).
- By x, the median is
(6,5). It becomes the root, splitting atx = 6. Left:(1,6) (3,2) (4,9). Right:(8,1) (9,7) (11,4). - The left half by y is
(3,2) (1,6) (4,9). Its median(1,6)splits aty = 6, with(3,2)below and(4,9)above. - The right half by y is
(8,1) (11,4) (9,7). Its median(11,4)splits aty = 4, with(8,1)below and(9,7)above.
Now ask for every point with 7 <= x <= 10 and 0 <= y <= 5. At the root, (6,5) is not inside. The query's left edge x = 7 is to the right of the split x = 6, so nothing in the left subtree can be inside: three points are skipped without being looked at. On the right, (11,4) is not inside; the query spans y = 4, so both children are visited, and (8,1) is the answer. Four nodes visited out of seven; with a million points the same pruning skips almost everything.
A k-d tree in one array
The tree does not need nodes and pointers. Put the median of the range [lo, hi) at position mid = lo + (hi - lo) / 2, the smaller points before it and the larger after it, and recurse on [lo, mid) and [mid + 1, hi). The array itself is the tree: the node for a range is the element in its middle. std::nth_element does exactly this step: it puts the element that belongs at mid there, nothing larger before it and nothing smaller after it, in linear time on average (the complexity the C++ standard specifies). Each level of the tree does O(n) work in total and there are about log2(n) levels, so building takes O(n log n) expected time and O(n) memory.
The comparator breaks ties on the splitting coordinate with the other coordinate. Without that, when several points share the median x, the standard leaves it unspecified which of them lands in the middle, and libstdc++ and libc++ may choose differently, so the tree's shape and every count below would depend on the compiler.
Ties still exist across the split: points with the same x as the splitting point can sit on either side. So the left subtree holds coordinates <= the split and the right subtree >= it, and the query must go left when the box reaches down to the split and right when it reaches up to it, both inclusive.
#include <algorithm>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
struct Box { long long xlo, ylo, xhi, yhi; }; // closed rectangle
bool inside(Point p, Box b) {
return b.xlo <= p.x && p.x <= b.xhi && b.ylo <= p.y && p.y <= b.yhi;
}
// Order by one axis, ties broken by the other, so every split is unique.
bool lessOn(int axis, Point a, Point b) {
return axis == 0 ? (a.x != b.x ? a.x < b.x : a.y < b.y)
: (a.y != b.y ? a.y < b.y : a.x < b.x);
}
// Implicit k-d tree: the node for [lo, hi) is pts_[mid]; its children are
// [lo, mid) and [mid + 1, hi). Even depths split on x, odd depths on y.
class KdTree {
public:
explicit KdTree(std::vector<Point> p) : pts_(std::move(p)) {
build(0, pts_.size(), 0);
}
// Appends every point inside b to out; returns the number of nodes visited.
long range(Box b, std::vector<Point>& out) const {
return range(0, pts_.size(), 0, b, out);
}
private:
void build(size_t lo, size_t hi, int axis) {
if (hi - lo <= 1) return;
size_t mid = lo + (hi - lo) / 2;
std::nth_element(pts_.begin() + lo, pts_.begin() + mid, pts_.begin() + hi,
[axis](Point a, Point b) { return lessOn(axis, a, b); });
build(lo, mid, 1 - axis);
build(mid + 1, hi, 1 - axis);
}
long range(size_t lo, size_t hi, int axis, Box b, std::vector<Point>& out) const {
if (lo >= hi) return 0;
size_t mid = lo + (hi - lo) / 2;
Point p = pts_[mid];
if (inside(p, b)) out.push_back(p);
long long split = axis == 0 ? p.x : p.y;
long long boxLo = axis == 0 ? b.xlo : b.ylo;
long long boxHi = axis == 0 ? b.xhi : b.yhi;
long visited = 1;
// The left part holds coordinates <= split, the right part >= split.
if (boxLo <= split) visited += range(lo, mid, 1 - axis, b, out);
if (boxHi >= split) visited += range(mid + 1, hi, 1 - axis, b, out);
return visited;
}
std::vector<Point> pts_;
};
int main() {
Lcg rng{42};
const Box query{400'000, 400'000, 410'000, 410'000}; // 0.01% of the area
std::println("{:>8} {:>6} {:>14} {:>12}", "n", "found", "nodes visited",
"brute found");
for (int n : {1'000, 10'000, 100'000, 1'000'000}) {
std::vector<Point> pts;
for (int i = 0; i < n; ++i)
pts.push_back({rng.range(0, 1'000'000), rng.range(0, 1'000'000)});
long brute = std::ranges::count_if(pts,
[&](Point p) { return inside(p, query); });
KdTree tree(pts);
std::vector<Point> found;
long visited = tree.range(query, found);
std::println("{:>8} {:>6} {:>14} {:>12}", n, found.size(), visited, brute);
}
}Press Run to compile and run this example.
n found nodes visited brute found
1000 0 10 0
10000 1 15 1
100000 20 54 20
1000000 107 194 107A thousandfold increase in n raised the work from 10 nodes to 194, and most of those 194 are the 107 points that had to be reported anyway. The guarantee behind this, for a balanced 2D k-d tree, is a worst case of O(sqrt(n) + k) nodes for a rectangle that contains k points (Lee and Wong, 1977). The sqrt(n) term counts cells that the rectangle's border cuts through: a vertical line can cross at most O(sqrt(n)) of the tree's cells, because every second level splits the other way and stops the line from passing through both children.
Seeing the partition
The program below builds the same tree for 40 points and draws each splitting line across the cell it divides, vertical for x-splits and horizontal for y-splits. The blue rectangle is a range query; the points it finds are red. The drawing helper flips the y-axis, because SVG's y grows downward.
#include <algorithm>
#include <format>
#include <print>
#include <string>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
class Svg {
public:
Svg(double width, double height) : w_(width), h_(height) {}
void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
a.x, fy(a.y), b.x, fy(b.y), color, width);
}
void dot(Point p, const char* color = "#1f2937", double r = 3) {
add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
}
void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
const char* fill = "none") {
std::string list;
for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
}
void text(Point p, const std::string& s, double size = 8) {
add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
}
void print() const {
std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
std::print("{}", body_);
std::println("</svg>");
}
private:
double fy(double y) const { return h_ - y; } // SVG's y axis points down
template <class... Args>
void add(std::format_string<Args...> fmt, Args&&... args) {
body_ += std::format(fmt, std::forward<Args>(args)...);
body_ += '\n';
}
double w_, h_;
std::string body_;
};
struct Box { long long xlo, ylo, xhi, yhi; };
bool inside(Point p, Box b) {
return b.xlo <= p.x && p.x <= b.xhi && b.ylo <= p.y && p.y <= b.yhi;
}
bool lessOn(int axis, Point a, Point b) {
return axis == 0 ? (a.x != b.x ? a.x < b.x : a.y < b.y)
: (a.y != b.y ? a.y < b.y : a.x < b.x);
}
// Builds the implicit tree and draws each split line across its cell.
void buildAndDraw(std::vector<Point>& pts, size_t lo, size_t hi, int axis, Box cell,
Svg& svg) {
if (lo >= hi) return;
size_t mid = lo + (hi - lo) / 2;
std::nth_element(pts.begin() + lo, pts.begin() + mid, pts.begin() + hi,
[axis](Point a, Point b) { return lessOn(axis, a, b); });
Point p = pts[mid];
Box left = cell, right = cell;
if (axis == 0) {
svg.line({p.x, cell.ylo}, {p.x, cell.yhi}, "#1f2937", 1);
left.xhi = right.xlo = p.x;
} else {
svg.line({cell.xlo, p.y}, {cell.xhi, p.y}, "#1f2937", 1);
left.yhi = right.ylo = p.y;
}
buildAndDraw(pts, lo, mid, 1 - axis, left, svg);
buildAndDraw(pts, mid + 1, hi, 1 - axis, right, svg);
}
int main() {
Lcg rng{11};
std::vector<Point> pts;
for (int i = 0; i < 40; ++i) pts.push_back({rng.range(5, 295), rng.range(5, 195)});
Svg svg(300, 200);
const Box world{0, 0, 300, 200}, query{150, 60, 250, 140};
svg.polygon({{query.xlo, query.ylo}, {query.xhi, query.ylo},
{query.xhi, query.yhi}, {query.xlo, query.yhi}},
"#2563eb", "rgba(37,99,235,0.15)");
buildAndDraw(pts, 0, pts.size(), 0, world, svg);
for (Point p : pts) svg.dot(p, inside(p, query) ? "#dc2626" : "#1f2937");
svg.print();
}Press Run to compile and run this example.
Every point sits on a line, because every point is a splitting point. Notice how the cells adapt to the data: where points are dense the cells are small, where they are sparse the cells are large. The query only has to look into cells that the blue rectangle overlaps, and the long first lines let it discard large parts of the picture at once.
Nearest neighbor: search near, then prune far
The nearest-neighbor search is a depth-first walk with a running answer best, the smallest squared distance found so far. At each node it measures the node's point, then descends first into the child on the same side of the split as the query q, because that is where the nearest point most likely is. When it comes back, it decides whether the other child is worth visiting. Every point on the other side of the split is at least |gap| away from q, where gap is the distance from q to the splitting line, so if gap * gap >= best, the whole subtree is skipped.
That test only uses one line. A stronger one keeps track of the rectangle (the cell) that each subtree covers and skips the subtree when the distance from q to the nearest point of that rectangle is already at least best. The cell test is what Friedman, Bentley and Finkel used in their 1977 paper. The program runs both against brute force, on uniform points and on points packed into 20 small clusters, with query points spread uniformly over the square.
#include <algorithm>
#include <limits>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
constexpr long long kMaxCoord = 1'000'000'000;
static_assert(8 * kMaxCoord * kMaxCoord < std::numeric_limits<long long>::max(),
"squared distances between points within kMaxCoord must fit");
long long dist2(Point a, Point b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
bool lessOn(int axis, Point a, Point b) {
return axis == 0 ? (a.x != b.x ? a.x < b.x : a.y < b.y)
: (a.y != b.y ? a.y < b.y : a.x < b.x);
}
struct Box { long long xlo, ylo, xhi, yhi; };
// Squared distance from q to the closest point of the rectangle b.
long long boxDist2(Point q, Box b) {
long long dx = std::max({b.xlo - q.x, 0LL, q.x - b.xhi});
long long dy = std::max({b.ylo - q.y, 0LL, q.y - b.yhi});
return dx * dx + dy * dy;
}
class KdTree {
public:
explicit KdTree(std::vector<Point> p) : pts_(std::move(p)) {
build(0, pts_.size(), 0);
}
// Squared distance from q to its nearest point. Adds visited nodes to count.
long long nearest(Point q, bool pruneByCell, long& count) const {
long long best = std::numeric_limits<long long>::max();
Box world{-kMaxCoord, -kMaxCoord, kMaxCoord, kMaxCoord};
nearest(0, pts_.size(), 0, world, q, pruneByCell, best, count);
return best;
}
private:
void build(size_t lo, size_t hi, int axis) {
if (hi - lo <= 1) return;
size_t mid = lo + (hi - lo) / 2;
std::nth_element(pts_.begin() + lo, pts_.begin() + mid, pts_.begin() + hi,
[axis](Point a, Point b) { return lessOn(axis, a, b); });
build(lo, mid, 1 - axis);
build(mid + 1, hi, 1 - axis);
}
void nearest(size_t lo, size_t hi, int axis, Box cell, Point q, bool pruneByCell,
long long& best, long& count) const {
if (lo >= hi) return;
if (pruneByCell && boxDist2(q, cell) >= best) return; // nothing here can win
size_t mid = lo + (hi - lo) / 2;
Point p = pts_[mid];
++count;
best = std::min(best, dist2(p, q));
long long gap = axis == 0 ? q.x - p.x : q.y - p.y; // signed, to the split
Box left = cell, right = cell;
if (axis == 0) left.xhi = right.xlo = p.x;
else left.yhi = right.ylo = p.y;
bool qLeft = gap < 0;
// Search the side that contains q first: it usually shrinks best quickly.
if (qLeft) nearest(lo, mid, 1 - axis, left, q, pruneByCell, best, count);
else nearest(mid + 1, hi, 1 - axis, right, q, pruneByCell, best, count);
// Every point on the other side is at least |gap| away from q.
if (gap * gap >= best) return;
if (qLeft) nearest(mid + 1, hi, 1 - axis, right, q, pruneByCell, best, count);
else nearest(lo, mid, 1 - axis, left, q, pruneByCell, best, count);
}
std::vector<Point> pts_;
};
// Uniform points in a square, or the same number in 20 small clusters.
std::vector<Point> makePoints(int n, bool clustered, Lcg& rng) {
std::vector<Point> centers, pts;
for (int i = 0; i < 20; ++i)
centers.push_back({rng.range(0, 1'000'000), rng.range(0, 1'000'000)});
for (int i = 0; i < n; ++i) {
Point c = clustered ? centers[rng.range(0, 19)] : Point{500'000, 500'000};
int r = clustered ? 2'000 : 500'000;
pts.push_back({c.x + rng.range(-r, r), c.y + rng.range(-r, r)});
}
return pts;
}
int main() {
Lcg rng{7};
const int queries = 200;
std::println("{:>9} {:>7} {:>12} {:>12} {:>6}", "data", "n", "split test",
"cell test", "wrong");
for (bool clustered : {false, true}) {
for (int n : {1'000, 10'000, 100'000}) {
std::vector<Point> pts = makePoints(n, clustered, rng);
KdTree tree(pts);
long bySplit = 0, byCell = 0, wrong = 0;
for (int i = 0; i < queries; ++i) {
Point q{rng.range(0, 1'000'000), rng.range(0, 1'000'000)};
long long brute = std::numeric_limits<long long>::max();
for (Point p : pts) brute = std::min(brute, dist2(p, q));
wrong += tree.nearest(q, false, bySplit) != brute;
wrong += tree.nearest(q, true, byCell) != brute;
}
std::println("{:>9} {:>7} {:>12.1f} {:>12.1f} {:>6}",
clustered ? "clustered" : "uniform", n,
double(bySplit) / queries, double(byCell) / queries, wrong);
}
}
}Press Run to compile and run this example.
data n split test cell test wrong
uniform 1000 15.1 14.9 0
uniform 10000 20.2 19.9 0
uniform 100000 23.4 23.1 0
clustered 1000 82.8 50.3 0
clustered 10000 642.2 172.6 0
clustered 100000 5781.3 544.3 0The columns are average nodes visited per query. On uniform data both tests behave like O(log n): a hundredfold more points cost about eight more nodes. That is the "logarithmic expected time" of the 1977 paper's title, and it holds for well-spread data, not for every input.
The clustered rows are the surprise, and they are worth understanding. A query point that lands in empty space is far from every cluster, so the circle around it with radius sqrt(best) is large, and its edge passes through the nearest cluster, where thousands of points are almost exactly as far away as the true nearest one. Every cell that the circle's edge crosses might contain a closer point and has to be opened. With the split test, the line distance is a poor bound and the cost grows almost linearly; with the cell test it grows roughly like sqrt(n). In both cases the answer is still exact. Only the cost changes.
Two facts follow. First, prune with the cell, not only the line; it costs four extra numbers per call. Second, the worst case of an exact nearest-neighbor query in a k-d tree is linear: put all points on a circle around the query and every point is a near-tie. Queries that come from the same distribution as the data, such as "the nearest driver to a customer in a city", are the friendly case; the Try it section lets you measure it.
Testing against brute force, with ties on purpose
The program above compared its answers with brute force, but on random coordinates up to a million, where two points almost never share an x. The bugs in spatial indexes live in the ties: a point exactly on a splitting line, a query box whose edge is exactly at a split, many identical points. So the test below uses tiny coordinate ranges, where ties are everywhere, plus whole sets of identical points and of collinear points. It compares range counts and nearest distances with a loop over all points. For nearest neighbor it compares the distance, not the point, because when two points are equally near either answer is correct.
#include <algorithm>
#include <limits>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
struct Box { long long xlo, ylo, xhi, yhi; };
bool inside(Point p, Box b) {
return b.xlo <= p.x && p.x <= b.xhi && b.ylo <= p.y && p.y <= b.yhi;
}
long long dist2(Point a, Point b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
bool lessOn(int axis, Point a, Point b) {
return axis == 0 ? (a.x != b.x ? a.x < b.x : a.y < b.y)
: (a.y != b.y ? a.y < b.y : a.x < b.x);
}
class KdTree {
public:
explicit KdTree(std::vector<Point> p) : pts_(std::move(p)) {
build(0, pts_.size(), 0);
}
long countIn(Box b) const { return countIn(0, pts_.size(), 0, b); }
long long nearest(Point q) const {
long long best = std::numeric_limits<long long>::max();
nearest(0, pts_.size(), 0, q, best);
return best;
}
private:
void build(size_t lo, size_t hi, int axis) {
if (hi - lo <= 1) return;
size_t mid = lo + (hi - lo) / 2;
std::nth_element(pts_.begin() + lo, pts_.begin() + mid, pts_.begin() + hi,
[axis](Point a, Point b) { return lessOn(axis, a, b); });
build(lo, mid, 1 - axis);
build(mid + 1, hi, 1 - axis);
}
long countIn(size_t lo, size_t hi, int axis, Box b) const {
if (lo >= hi) return 0;
size_t mid = lo + (hi - lo) / 2;
Point p = pts_[mid];
long long split = axis == 0 ? p.x : p.y;
long n = inside(p, b);
if ((axis == 0 ? b.xlo : b.ylo) <= split) n += countIn(lo, mid, 1 - axis, b);
if ((axis == 0 ? b.xhi : b.yhi) >= split) n += countIn(mid + 1, hi, 1 - axis, b);
return n;
}
void nearest(size_t lo, size_t hi, int axis, Point q, long long& best) const {
if (lo >= hi) return;
size_t mid = lo + (hi - lo) / 2;
Point p = pts_[mid];
best = std::min(best, dist2(p, q));
long long gap = axis == 0 ? q.x - p.x : q.y - p.y;
if (gap < 0) nearest(lo, mid, 1 - axis, q, best);
else nearest(mid + 1, hi, 1 - axis, q, best);
if (gap * gap >= best) return;
if (gap < 0) nearest(mid + 1, hi, 1 - axis, q, best);
else nearest(lo, mid, 1 - axis, q, best);
}
std::vector<Point> pts_;
};
int main() {
Lcg rng{2024};
long cases = 0, failed = 0;
for (int trial = 0; trial < 3'000; ++trial) {
int n = rng.range(1, 60);
// Sides of 3 and 10 make equal coordinates common; 1,000 makes them rare.
int side = trial % 3 == 0 ? 3 : trial % 3 == 1 ? 10 : 1'000;
std::vector<Point> pts;
for (int i = 0; i < n; ++i)
pts.push_back({rng.range(0, side), rng.range(0, side)});
if (trial % 10 == 0) pts.assign(n, pts[0]); // all identical
if (trial % 10 == 1) for (Point& p : pts) p.y = p.x; // all on one line
KdTree tree(pts);
for (int q = 0; q < 5; ++q) {
int x0 = rng.range(-1, side + 1), y0 = rng.range(-1, side + 1);
Box b{x0, y0, x0 + rng.range(0, side), y0 + rng.range(0, side)};
Point c{rng.range(-1, side + 1), rng.range(-1, side + 1)};
long inBox = std::ranges::count_if(pts,
[&](Point p) { return inside(p, b); });
long long near = std::numeric_limits<long long>::max();
for (Point p : pts) near = std::min(near, dist2(p, c));
++cases;
failed += tree.countIn(b) != inBox || tree.nearest(c) != near;
}
}
std::println("{} queries checked against brute force, {} failed", cases, failed);
}Press Run to compile and run this example.
15000 queries checked against brute force, 0 failedA passing test is only convincing if it can fail. Change the first <= split in countIn to < split, a plausible slip, and the test reports 2467 failures. Every one of them comes from the trials with sides 3 and 10; the trials with side 1,000 catch none. A test on "realistic" random data would have passed the broken tree.
Point-region quadtrees: split space, not data
A k-d tree splits at data points. A point-region (PR) quadtree splits space: the root is a square, and a square that holds too many points is cut into four equal quadrants, whatever the points look like. (The quad tree of Finkel and Bentley, 1974, split at data points, like a k-d tree with four children; the point-region variant always halves the square.)
Fixed splits make updates easy: an insert walks down to the leaf whose square contains the point and splits that leaf if it overflows, with no medians and no rebalancing. Two numbers control it: the bucket capacity, how many points a leaf holds before it splits, and the maximum depth.
The maximum depth is not optional, because the depth depends on how close points are, not on how many there are: for a quadtree that splits until each leaf holds one point, the depth is at most log2(s / c) + 3/2, where s is the side of the root square and c the smallest distance between two points (de Berg et al., Lemma 14.1). Two points one unit apart in a square of side 2^20 force about 20 levels; two identical points force infinitely many. The program builds three trees (2,000 uniform points; 1,000 uniform points plus five points one unit apart; 2,000 copies of one point) and checks a range query on each against brute force.
#include <algorithm>
#include <array>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
struct Box { long long xlo, ylo, xhi, yhi; };
bool inside(Point p, Box b) {
return b.xlo <= p.x && p.x <= b.xhi && b.ylo <= p.y && p.y <= b.yhi;
}
// Point-region quadtree over the square [0, size) x [0, size).
// A leaf splits into four equal quadrants when it holds more than kCapacity
// points, unless it is already at kMaxDepth: then it simply keeps them all.
class QuadTree {
public:
static constexpr int kCapacity = 4, kMaxDepth = 16;
explicit QuadTree(long long size) { nodes_.push_back({0, 0, size, 0, {}, {}}); }
void insert(Point p) {
int i = 0;
while (nodes_[i].child[0] != 0) i = nodes_[i].child[quadrant(nodes_[i], p)];
nodes_[i].pts.push_back(p);
Node& leaf = nodes_[i];
if (int(leaf.pts.size()) > kCapacity && leaf.depth < kMaxDepth) split(i);
}
long countIn(Box b, long& visited) const { return countIn(0, b, visited); }
int depth() const {
int d = 0;
for (const Node& n : nodes_) d = std::max(d, n.depth);
return d;
}
size_t largestLeaf() const {
size_t m = 0;
for (const Node& n : nodes_) m = std::max(m, n.pts.size());
return m;
}
size_t nodeCount() const { return nodes_.size(); }
private:
struct Node {
long long x0, y0, size;
int depth;
std::array<int, 4> child; // indices into nodes_; all 0 for a leaf
std::vector<Point> pts; // only leaves hold points
};
static int quadrant(const Node& n, Point p) {
long long h = n.size / 2;
return (p.x >= n.x0 + h ? 1 : 0) + (p.y >= n.y0 + h ? 2 : 0);
}
void split(int i) {
long long h = nodes_[i].size / 2;
for (int q = 0; q < 4; ++q) {
Node c{nodes_[i].x0 + (q & 1) * h, nodes_[i].y0 + (q >> 1) * h, h,
nodes_[i].depth + 1, {}, {}};
nodes_.push_back(c); // may reallocate: hold indices, not references
nodes_[i].child[q] = int(nodes_.size()) - 1;
}
std::vector<Point> moved = std::move(nodes_[i].pts);
nodes_[i].pts.clear();
for (Point p : moved) insert(p);
}
long countIn(int i, Box b, long& visited) const {
const Node& n = nodes_[i];
++visited;
bool misses = b.xhi < n.x0 || b.xlo >= n.x0 + n.size ||
b.yhi < n.y0 || b.ylo >= n.y0 + n.size;
if (misses) return 0;
long found = std::ranges::count_if(n.pts, [&](Point p) { return inside(p, b); });
if (n.child[0] != 0)
for (int c : n.child) found += countIn(c, b, visited);
return found;
}
std::vector<Node> nodes_;
};
void report(const char* name, const std::vector<Point>& pts) {
QuadTree tree(1 << 20);
for (Point p : pts) tree.insert(p);
Box q{400'000, 400'000, 420'000, 420'000};
long visited = 0;
long found = tree.countIn(q, visited);
long brute = std::ranges::count_if(pts, [&](Point p) { return inside(p, q); });
std::println("{:<14} nodes {:>4} depth {:>2} largest leaf {:>4} found {:>4}/{:<4}"
" visited {}", name, tree.nodeCount(), tree.depth(), tree.largestLeaf(),
found, brute, visited);
}
int main() {
Lcg rng{5};
std::vector<Point> uniform, close, same;
for (int i = 0; i < 2'000; ++i)
uniform.push_back({rng.range(0, (1 << 20) - 1), rng.range(0, (1 << 20) - 1)});
close.assign(uniform.begin(), uniform.begin() + 1'000);
for (int i = 0; i < 5; ++i) close.push_back({410'000 + i, 410'000}); // 1 unit apart
same.assign(2'000, Point{410'000, 410'000});
report("uniform", uniform);
report("close points", close);
report("all identical", same);
}Press Run to compile and run this example.
uniform nodes 1465 depth 7 largest leaf 4 found 1/1 visited 21
close points nodes 773 depth 16 largest leaf 5 found 5/5 visited 65
all identical nodes 65 depth 16 largest leaf 2000 found 2000/2000 visited 65Uniform data gives depth 7. Five points within four units of each other drive a chain of splits down to the cap; 2,000 identical points would split forever without it, because every split sends all of them into the same quadrant. With the cap, the tree ends in one large leaf that a query scans linearly, which is the honest cost for data that really is all in one place. The node counts show the other price of splitting space: the identical points alone create 65 nodes, four per level, almost all of them empty. A k-d tree over n points has exactly n nodes.