Chan Alston | Full stack developer

For every one of the ~7,000 detections in an image, I needed a statistic computed over that detection’s ~120 nearest neighbours. Doing it naively means measuring every point against every other point, roughly 49 million distance evaluations per image, with several images running concurrently. So, obviously: build a spatial index.

I built a uniform grid, and it worked. Then I tried tuning K for a different data sets that have different feature to detect. The result is the same, nothing changes. It feels like somewhere past a few hundred, K simply stopped doing anything, and no error told me.

The problem

Okay, say you are handed a set of points in a 2D image, say ~7,000 detections. The distribution is roughly uniform, with each point having a feature attached to it. You need to compare each point’s feature to the features of its K nearest neighbours to compute a local statistic.

Seems simple enough, right? Just measure the distance for each point, and compute the statistic. However, that would be too slow because each point has to compare itself to every other point, which is O(n²) in time complexity. Instead, we can use a spatial data structure to speed up the nearest neighbour search.

The query set is the point set is an important point. Building the index accounts for about O(n) while the query is about O(K) which is much smaller than brute force with O(n²) of query.

Spatial Grid to the Rescue

A spatial grid is a simple and effective way to partition the 2D space into cells, allowing us to quickly find nearby points. By dividing the image into a grid of cells, we can limit our search for nearest neighbours to only those points that are in the same cell or nearby cells, widening the search ring by ring around the query’s cell until K neighbours have been found.

You can play with the parameters of the spatial grid and see how it affects the nearest neighbour search in the widget below. Adjust the cell size and K value to see how they impact the number of neighbours found and the performance of the algorithm.

The grid answers “who is nearby?” without visiting every point

144 detections on an X–Y plane, bucketed into an N×N grid of cells. Hover any point to query it.

K — neighbours requested10

Grid size — 16.7 units/cell · 4.0 pts/cell6

Neighbours found10 / 10

Rings used1 / 4

Candidates examined38

K selected candidates examined untouched

The search widens ring by ring around the query's cell and stops as soon as K neighbours are found — or gives up after ring 4. Push K high enough and watch the search come back short: the geometry of the index, not your parameter, decides what is reachable.

There is a catch, though. The cell size is a critical parameter that determines how many points are in each cell.

If the cell size is too small, we may not find enough neighbours (even with ring expansion), which silently limits the maximum K value we can use.
If the cell size is too large, we may end up checking too many points, which can slow down the search.

This can be an issue: if the cell size is not chosen correctly, changing the K parameter may have no effect on the results. No crash, no NaN, no error bar screaming at you — just the same number, quietly wrong.

Why not a KD-tree?

Using a KD-tree is another common approach for nearest neighbour search, but it isn’t the right fit here:

The distribution is roughly uniform and bounded by the image rectangle. A KD-tree’s adaptive splits pay off on clustered or unbounded data; with nothing to adapt to, its extra structure is wasted.
Query set == point set, so we build the index once and run n queries either way; The O(n log n) KD-tree build and the O(n) grid build are both amortized over those n queries. The grid just has the cheaper build.
The grid’s query is expected O(K) distance evaluations that walk a contiguous CSR array, while a KD-tree query backtracks through pointer linked nodes. For large K that pointer chasing and poor cache locality dominate, whereas the grid’s flat, sequential access stays fast.

KD-trees aren’t bad, they’re just not the right tool for this specific problem.

What about Quadtree?

Quadtree is another spatial data structure that can be used for nearest neighbour search. However, it has similar issues to KD-trees in this case. The uniform distribution of points and the bounded nature of the image leave the quadtree’s adaptive subdivision with nothing to adapt to, all these extra gymnastics are wasted.

If the points were clustered or had varying density, a quadtree might be a better choice. But in this case, the spatial grid is the most appropriate data structure for the problem at hand.

The failure mode: my first guess capped K

My first answer to “so, what is the right cell size?” was the obvious one - the kind of answer that feels right in the moment and bites you later. The typical gap between points, the pitch, can be estimated from the bounding box alone (the derivation is coming), so make each cell two pitches wide. It looks fine where there are roughly four points per cell, and most queries finish within ring 1.

But a ring-limited search has a hard ceiling that you can write down. A search that gives up after maxRing rings can reach at most:

K_{\text{max}} = (2 \cdot \text{maxRing} + 1)^2 \cdot \left(\frac{\text{cell}}{\text{pitch}}\right)^2

With maxRing = 4 and cell = 2 · pitch: K_max = 81 * 4 = 324. Now replay what happens:

Ask for K ≤ 324: everything works, and the results look reasonable.
Ask for K = 500: the rings exhaust themselves, the loop exits at maxRing, and the search returns the 324 neighbours it managed to gather. No error, no warning.
The statistics are computed from 324 neighbours whether you asked for 400 or 4,000.

The parameter did not fail. It was capped by the geometry of the index, and the number I typed only agreed with reality by coincidence.

The cell size must be computed backwards from K.

The solution: compute the cell size based on K

Everything we need is already in the data. Start with the bounding box, the largest and smallest coordinates:

w = x_{\text{max}} - x_{\text{min}}, \qquad h = y_{\text{max}} - y_{\text{min}}

Assuming the points are uniformly distributed across that area, every point owns roughly the same amount of real estate: a patch of $\frac{w \cdot h}{n}$ . Note what is not happening here; no distance between points is ever measured. We take the side length of that patch as an estimate of the typical gap between neighbouring points, and call it the pitch:

\text{pitch} = \sqrt{\frac{w \cdot h}{n}}

For detections laid out on a regular pattern, dots on a grid, LEDs on a panel, that patch is close to a literal square, and its side really is the centre to centre spacing, so the estimate lands almost exact. For noisier layouts the number drifts away from any true pairwise distance; treat it as a scale, not a measurement. Either way, it is precisely the quantity the cell size rule needs.

Now, how many points does a cell of side length $\text{cell}$ hold? A cell spans $\left(\frac{\text{cell}}{\text{pitch}}\right)^2$ pitch sized squares worth of area, so on average:

\text{pts per cell} = \left(\frac{\text{cell}}{\text{pitch}}\right)^2

A search that expands out to ring $r$ scans a $(2r+1)^2$ block of cells, so the number of points reachable at ring $r$ is:

\text{capacity}(r) = (2r+1)^2 \cdot \left(\frac{\text{cell}}{\text{pitch}}\right)^2

Here is the requirement that drives everything: ring 1 alone - the $3 \times 3$ block of 9 cells, should comfortably hold $K$ . We target $2K$ rather than $K$ for slack against local density variation:

9 \cdot \left(\frac{\text{cell}}{\text{pitch}}\right)^2 = 2K

Solve for the ratio, then scale it by the pitch to get the actual cell dimension:

\frac{\text{cell}}{\text{pitch}} = \sqrt{\frac{2K}{9}} \qquad \Longrightarrow \qquad \boxed{\;\text{cell} = \text{pitch} \cdot \sqrt{\frac{2K}{9}}\;}

Sanity check with $K = 120$ : $\sqrt{2 \cdot 120 / 9} \approx 5.16$ , so the cell is about five pitches wide, and ring-1 capacity is $9 \cdot 5.16^2 \approx 240 = 2K$ . ✓

When I wrote it out this way the fix felt obvious, but getting here took longer than I’d like to admit. The entire derivation is just algebra — but the hard part was realizing that the cell size was my problem to solve, not a free parameter I could set once and forget.

One clarification worth pinning down: $2K$ is the size of the candidate pool that ring 1 is sized to hold, not the amount returned. Every candidate in the pool competes by distance, and the search keeps exactly the best $K$ of them (or fewer, plus the not judgeable sentinel). The doubling exists purely so that this cut usually happens while still inside ring 1.

Why target 2K and not just K?

Two things to take note here. The first is plain slack against local density variation: real detections clump, and when a stretch of the image runs twice as dense as the global estimate, those regions still turn up plenty of neighbours inside ring 1.

The second is subtler. The scanned region is a square of cells, but the true $K$ -th-neighbour boundary is a disc, but we are using squares here. This mismatch is the source of a blind spot: a nearer point can sit just outside the block, never scanned, while a farther insider takes its slot.

For a query sitting near the middle of its cell, ring-1’s half width is $1.5 \cdot \text{cell} = 1.5\sqrt{2K/9} \cdot \text{pitch}$ , while the mean field $K$ -th distance is $\sqrt{K/\pi} \cdot \text{pitch}$ , a ratio of $\frac{\sqrt{2\pi}}{2} \approx 1.25$ , independent of K. That margin shrinks as the query drifts toward a corner of its cell, bottoming out near $0.84$ .

This solves the issue completely. Under the fixed rule, ring-1 capacity was some constant the cell size happened to imply, and $K$ was silently capped by it. Under this rule, ring-1 capacity is $2K$ for any $K$ , so the search is never ring-limited and $K$ always means what it says.

Two practical floors are applied on top:

$\frac{\text{cell}}{\text{pitch}} \geq 2.0$ , binds below $K = 18$ , preventing tiny values of $K$ from producing absurdly small cells and a huge grid array.
$\text{cell} \geq 4\,\text{px}$ , guards against degenerate bounding boxes and sub-pixel pitches.

And yes - the uniformity assumption is doing work here. That is precisely what the factor of two is insurance against.

The numbers

Operation counts, not wall clock; At these sizes wall clock is dominated by whatever else the pipeline is doing.

	Brute force	CSR grid
Build	-	`O(n)`, two linear passes
Per query	`n` distance evaluations	≈ `2K` candidates examined
`n = 7,000`, `K = 120`	~7,000	~240 → ~29x fewer
Total, `n = 7,000`	~49M candidate touches	~1.7M candidate touches

One non obvious consequence of the sizing rule: expected points per cell is $2K/9$ , so the offset table holds only $n / (2K/9)$ cells, about $n/27$ at $K = 120$ . The index is smaller than the point array, and it shrinks further as K grows.

Compressed Sparse Row (CSR)

To store the grid efficiently, we borrow the Compressed Sparse Row idea from sparse matrix formats: all points go into one flat items[] array, grouped contiguously by cell, with a start[] offset table holding one entry per cell (plus one). An empty cell costs a single integer, not a list object. The construction is a counting sort: tally each point into start[b + 1], prefix sum the tallies into offsets, then scatter points through a disposable clone of start. Two linear passes over the points, zero allocation per point.

The widget below demonstrates exactly that construction. Step through it and watch the arrays fill.

How CSR gets built

Two passes over the points and one sweep of additions — zero allocation per point. Click step, or press play.

cellOf

start

items

ready

Seventeen detections scattered unevenly across a 4×4 grid — some cells empty, some crowded. Step through the three passes that turn them into flat arrays.

1/69

To answer a query, locate the point’s cell and collect candidates from it, then widen the search ring by ring until enough neighbours are found or the maximum ring limit is reached. Each ring restarts collection instead of extending the previous one, deliberate, since the common case exits at ring 1 and rescanning beats bookkeeping. Distances stay squared throughout (identical ordering, no sqrt in the innermost loop), and selection inserts into a fixed K length sorted array rather than sorting all candidates.

Boundaries

Ring expansion is limited by the image boundaries. If a point is near the edge of the image, the search may not be able to expand fully in all directions. However, it can still expand in the directions that are available until enough neighbours are found or the maximum ring limit is reached.

If the search reaches the maximum ring limit without finding enough neighbours, a sentinel value is returned to indicate that the query point is not judgeable. This allows the algorithm to handle cases where there are not enough neighbours without producing garbage statistics.

Yes, it is only approximation

The scanned region is a square of cells, but the true K nearest set is a disc. The 25% margin shown above keeps that mismatch rare, it does not make it impossible. A hot spot can still push a true neighbour outside the block, swapping it for a farther insider. Exactness would require an “expand until the ring’s inscribed circle exceeds the current K th distance” test, which is more work than this application justifies. Occasional swaps are acceptable noise, especially next to the alternative, a sentinel that admits the answer could be garbage.

Wrap up

This approach computes K nearest neighbours over large point sets with predictable behaviour: the cell size derives from the data, ring capacity tracks K by construction, and the index declines to answer when it genuinely cannot.

When a knob does nothing, the problem isn’t the knob, it’s something else entirely, and usually it’s the infrastructure around it silently eating your input. Some parameters shouldn’t be exposed to the user at all. The cell size is one of them: it should be computed from K and the data, not set arbitrarily and left to rot.

My K parameter did nothing