Thanks to @David Eisenstadt for the idea of a 3d search structure. This is part of the best answer, though his strange metric is not needed.
The key is to look in detail at how nearest neighbor search works. I'll show this for quadrees. Kd-trees with k=3 are similar. Here is pseudocode:
# Let nearest_info be a record containing the current nearest neighbor (or nil
# if none yet) and the distance from point to that nearest neighbor.
def find_nearest_neighbor(node, target, nearest_info)
if node is leaf
update nearest_info using target and the points found in this leaf
else
for each subdivision S of node
if S contains any point P where dist(P,T) < nearest_info.distance,
find_neareast(S, target, nearest_info)
end
end
end
end
When this is done, nearest_info
contains the nearest neighbor and its distance.
The key is if S contains any point P where dist(P,T) < nearest_info.distance
. In a 3d space, of (x,y,r)
triples that describe circles, we have
def dist(P,T)
return sqrt( (P.x - T.x)^2 + (P.y - T.y)^2 ) - P.r - T.r
end
Here P is an arbitrary point in an octant of an octree cuboid. How to consider all points in the cuboid? Note all components of T are effectively fixed for a given search, so it's clearer if we write the target as a constant point (a, b, c)
:
def dist(P)
return sqrt( (P.x - a)^2 + (P.y - b)^2 ) - P.r
end
Where we have left out c = T.r
completely because it can be subtracted out of the minimum distance after the algorithm is complete. In other words, the radius of the target does not affect the result.
With this it is pretty easy to see that the P
we need to obtain minimum dist to the cuboid is Euclidean closest to the target with respect to x and y and with the max represented radius. This is very easy and quick to compute: a 2d point-rectangle distance and a 1d max
operation.
In hindsight all this is obvious, but it took a while to see it from the right point of view. Thanks for the ideas.