9 K-d Trees: Multi-dimensional Data Indexing

Notice that in the Cartesian plane a vertical line through a point P=(Px,Py)
has a peculiar characteristic: it is parallel to the y
axis, and all points on the line have the same value for their x
coordinate, Px
. Likewise, a horizontal line passing through P
is made of all the points in the plane for which y=Py
.

Figure 9.2 Partitioning points in the 2-D Cartesian space by cycling through the directions along which we split. For the first split (left) we choose point R
and draw a vertical line (parallel to y
axis, x
coordinate is constant) passing through it. We have such created 2 half-spaces, on the left (yellow) and right (blue) of this line, grouping points W
, P
, O
, Q
, U
on one side, and S
, T
on the other. Point R
is the pivot of this partitioning.
Next, we choose point W
in the yellow partition: this time, we draw a horizontal line (parallel to x
axis, y
coordinate is constant): it splits the yellow partition into 2 new partitions, one in the top-left area of the plane (red), containing P
, O
and U
, and one in the bottom-left area (green), with just Q
.
If we further split the red area at point P
, we need to use again a vertical line, as shown in the right-most part of the figure.
So, when we split the plane along a vertical line through P=(Px,Py)
, what we really mean is that we create two partitions, one for the points L
in the plane for which Lx<Px
, and one for those points R
for which Rx>Px
. And similarly, for horizontal lines, using the y
coordinates of the points.

Can this algorithm be generalized to higher dimensions? Yes, it naturally allows a generalization to higher dimensions, because we split on a single coordinate for each point, but we can do a round-robin through all the coordinates of a k
-dimensional space, and the i
-th level of the binary tree we build will correspond to a split along the (i mod k)-
th dimension.
Meaning that in a 2-D space, the root will split the plane along the x
axis, its children will split each of the two semi-planes along the y
axis, then their children again along the x
axis and so on and so forth. In a k
-dimensional space, with k > 2
, we start with the first coordinate at level 0
(the root), then move to the second at height 1
, the third at height 2
etc...

(Left)     We add the first pivot, that is going to be the root of the tree. The pivot creates an implicit division along the x
axis, partitioning the remaining points into left and right subsets. The region covered by a node is the union of the partitions its pivot creates, so the root covers the whole subset.
(Center) Each of the two sub-regions implied by the root is further partitioned, by selecting a point in each region as pivot. As shown by the highlighting above the horizontal axis, now each of these nodes at level 1 covers half of the space (while the root still covers the whole space).
(Right) Level 2 nodes are added, further partitioning the space. Notice how some regions are still only covered by nodes at level 1, because these intermediate nodes only have one child.

A k-d tree is a binary search tree whose elements are points taken from a k
-dimensional space, i.e. tuples with k
elements, whose coordinates can be compared (to simplify, let’s assume each coordinate’s values can be translated into real numbers). In addition to that, in a k-d tree, at each level i
we only compare the i
-th (modulo k
) coordinate of points, to decide which branch of the tree will be traversed.

Figure 9.3 shows the construction of an example tree for the unidimensional case. This is an edge case, because singletons (tuples with dimension 1) will result in always using the x
coordinate in points (i.e. the whole singleton).

Here, when we use the term “cover”, we mean that, given an entry X
in the search space (in 1-D, given a real number X
), and if we query the tree (as a binary search tree) for X
, then all nodes we traverse during the search for X
(which form a path from the root to a node N
) cover X
. In particular, the node on which search stops is the one covering the smallest region for X
.

(Left)     First we split vertically with pivot “Opal City”, which becomes the root of the k-d tree. Vertical splits are drawn in red, as well as edges corresponding to vertical splits. The background of each tree node is filled with the colors of the regions the node covers, so the root is half green and half yellow, because it covers both partitions created with its pivot.
(Right)    The right partition (created by the root) is further split, along a horizontal line passing through our second pivot, “Civic City”, so in the BST we add a left child to the root. This node corresponds to another split into a top and bottom sub-region. Horizontal splits are drawn in blue, as well as corresponding edges in the tree.

There is not, however, a direct match, and looking at regions is not enough to determine the closest point: for instance, if you look at B-8
in the picture, it’s not clear if Buffalo, Pittsburgh or Harrisburg is the closest city, and C-6
looks closer to Syracuse than Harrisburg, despite being “covered” by the latter and not the former.

1. All points in the tree have dimension k
.

3. If a node N
’s split coordinate index is j
, then N
’s children have split coordinate equal to (j+1) mod k
.
4. For each node N
, with split coordinate index j
, all nodes L
in its left subtree have a smaller value for N
’s split coordinate, L[j] < N[j]
, and all nodes R
on N
’s right subtree have a larger or equal value for N’s split coordinate, R[j] ≥ N[j]
.

So far in this section we consciously ignored a detail, a special characteristic we look for when it comes to binary search trees: whether the tree is balanced or not. As for the regular BST
, you might have already figured out that the order in which we insert elements in our k-d tree determines the shape of the tree. Even having finite areas is not a synonym of having small areas, nor good partitioning: if you look at the example in figure 9.3, we have carefully chosen points “by hand” in order to create a balanced tree, but it is easy to provide an example of a terrible insertion sequence that would create an imbalanced tree (just insert points starting from the top-left going in order towards the bottom-right corner, for instance).

class KdNode
  #type tuple(k)
  point
  #type KdNode
  left
  #type KdNode
  right
  #type integer
  Level
  function KdNode(point, left, right, level)
  
class KdTree
  #type KdNode   root
  #type integer
  k
  function KdTree(points=[])

A KdTree
just contains a root node, and a constructor method, taking an optional array of points as input. We’ll take a look at how a k-d tree can be constructed later, after we introduce the insertion method. For now, it will suffice to say that it will set the root to either a void entry (being it null
, a special instance of KdNode
, or whatever it’s more appropriate for the language used).
For the sake of convenience, let’s assume a tree also stores the value k
, the dimension of the space on which the tree is defined.

In section 9.2 we have implicitly described how search works on k-d trees. There is nothing fancy in this algorithm, it’s just a regular search on a binary search trees storing tuples, with the caveat that instead of comparing the whole tuple at each step, we only use one coordinate: at level i
, we compare the i
-th coordinate (or i mod k
, if i≥k
).

function getNodeKey(node)                                               #1
  return getPointKey(node.point, node.level)                            #2
 
function getPointKey(point, level)                                      #3
  j ← level % k                                                         #4
  return point[j]                                                       #5
 
function compare(point, node)                                           #6
  return sign(getPointKey(point, node.level) - getNodeKey(node))        #7
 
function splitDistance(point, node)                                     #8
  return abs(getPointKey(point, node.level) - getNodeKey(node))         #9
 

Figure 9.9 An example of an unsuccessful search on a k-d tree (2-D). The searched point, P
, would ideally lie in the region highlighted in red, which corresponds to the left subtree of node D
.
Listing 9.3, instead, shows the pseudocode for the search method. The pseudo-code for all these methods will assume that we are providing an internal version that takes a KdNode
as argument. The public API for KdTree
, instead, will provide adapter methods that will call the internal ones, for instance KdTree::search(target)
will call search(root, target)
.

function search(node, target)                                           #1
  if node == null then                                                  #2
    return null
  elsif node.point == target then                                               #3
    return node
  elsif compare(target, node) < 0 then                                  #4
    return search(node.left     , target)
  else                                                          #5
    return search(node.right, target)
 

Figure 9.10 An example of a successful search on a k-d tree (2-D). Here points P
and E
are coincident, and the highlighted region corresponds to the subtree rooted at node E
.

function insert(node, newPoint, level=0)                                        #1
  if node == null then                                                  #2
    return new KdNode(newPoint, null, null, level)
  elsif node.point == newPoint then                                     #3
    return node
  elsif compare(newPoint, node) < 0 then                                        #4
    node.left ← insert(node.left, newPoint, node.level + 1)            #5
    return node
  else                                                          #6
    node.right ← insert(node.right, newPoint, node.level + 1)          #7
    return node
 

Let’s follow the first example step by step: it starts with a call to insert(A, (-1.5, 2))
, where we don’t pass any value for level
, thus defaulting it to the right value for root (as defined in the function signature, this value is 0
).
A <> null
, so condition at line #2 won’t match; A.point <> (-1.5, 2)
, and also at line #3 the condition is false. When we get to line #4, instead, -1.5 < 0
, so compare
will return -1
and we traverse the left subtree and call insert(C, (-1.5, -2), 1)
.
The next few calls will proceed similarly (like we have seen for search
), and in turn we call: insert(D, (-1.5, -2), 2), insert(null, (-1.5, -2), 3)
. For the latter, condition at line #2 will be true, so we create a new node, KdNode((-1.5, -2), null, null, 3)
, and return it to the previous call in the stack trace: there, at line #5, we set D.left
to this new KdNode
we created, and then return D
.

insert(A, (2.5, -3))
  insert(B, (-1.5, 2), 1)
      insert(E, (-1.5, 2), 2)
     insert(null, (-1.5, 2), 3)
     return new KdNode((-1.5, 2), null, null, 3)
    return E
  return B
return A

Before moving on to the most advanced methods designed for k-d trees, let’s take a step back. In section 9.2.3 we have already described one key point of k-d trees: we need our tree to be balanced. In the previous sections we have indeed seen how search and insertion have running time O(h)
, proportional to the height of the tree, and thus having a balanced tree will mean that h = log(n)
and in turn that all these methods will run in logarithmic time.

function compare(point, node)                                           #1
  s ← sign(getPointKey(point, node.level) - getNodeKey(node))          #2
  if s == 0 then                                                                #3
    return node.level % 2 == 0 ? -1 : +1
  else
    return s

At date, there is not a solution to easily keep a k-d tree balanced while maintaining O(h)
running time for insertion.
Nevertheless, if we know the set of points to insert in advance (when we create the k-d tree), then we can find an optimal order of insertion that allows us to construct a balanced tree. If the tree is not changed after construction, or if the number of elements inserted/deleted is negligible in comparison to the size of the tree, then we can have a worst-case balanced k-d tree, and insert
and search
will run in worst-case logarithmic time.

function constructKdTree(points, level=0)                               #1
  if size(points) == 0 then                                             #2
    return null
  elsif size(points) == 1 then                                          #3
    return new KdNode(points[0], null, null, level)
  else
    (median, left, right) ← partition(points, level)                   #4
      leftTree ← constructKdTree(left, level + 1)                      #5
      rightTree ← constructKdTree(right, level + 1)                    #6
    return new KdNode(median, leftTree, rightTree, level)               #7

The key point in listing 9.6 is the call to the partition
method at line #4: we need to pass level
as an argument because it will tell us which coordinate we need to use to compare points. The result of this call will be a tuple with the median of the points
array and two new arrays with each, (n-1)/2
elements, if size(points) == n
.
Each point in left
will be “smaller” (with respect to the coordinate at index level % k
) than median
, and each point in right
will be larger than median
: we can therefore recursively construct both (balanced) subtrees using these two sets.

To understand how it works, let’s consider a call to constructKdTree([(0,5),(1,-1),(-1,6),(-0.5,0),(2,5),(2.5,3),(-1, 1),(-1.5,-2)])
.
The median for this set (WRT the first coordinate, i.e. the median of all the first values of the tuples) is either -0.5
or 0
: there is an even number of elements, so there are technically 2 medians - you can double check the values by sorting the array.
Say we choose -0.5
as the median, then we have
(median, left, right) ←
(-0.5,0), [(-1, 1),(-1.5,-2),(-1,6)], [(1,-1), (2.5,3),(2,5)(0,5)]
So, at line #5 we call constructKdTree([(-1, 1),(-1.5,-2),(-1,6)], 1)
to create root’s left subtree. This in turn will partition again the sub-array, but comparing the second coordinates of each tuple: the median of y
coordinates is 1, so we have:
(median, left, right) ←
(-1, 1), [(-1.5,-2)], [(-1,6)]

Finally line #4, where we call partition
: it’s possible to find a median in linear time, and we can as well partition an array of n
elements around a pivot with O(n)
swaps (or create two new arrays, also with O(n)
total assignments).

To complete our analysis, if we look at the extra memory needed by this method, it will of course require O(n)
memory for the tree. There is, however, more: if we don’t partition the array in-place, and create a copy for each left and right sub-array, then each call to partition would use O(n)
extra memory, and deriving a similar formula to the one for the running time, we could find out that also M(n) = O(n * log(n))
.

That’s because we would need only a constant amount of memory for the internals of the function, plus the O(n)
memory needed for the tree.

After search
and insert
, we could only continue with the third basic operation on a container, remove
. This is despite the fact that on a k-d tree, delete is not such a common operation, and some implementations don’t even offer this method: as we discussed in previous section, k-d trees are not self-balancing, and so they perform best when they are created with a static set of points, avoiding frequent insertion and removal.

Nevertheless, in any real-world application you’ll likely need to be able to update your dataset, so we are going to describe how to remove elements: figures 9.13 and 9.14 shows the remove
method in action on our example k-d tree, and the result of removing point D
from it.
Similarly to insert
and search
, this method is based on binary search tree removal. There are, however, two issues that makes k-d tree’s version sensibly more complicated, and they are both connected to how we find a replacement in the tree for the node we are going to remove.

Figure 9.14 The k-d tree resulting after deleting point D
.
  A.        The node we need to remove is a leaf: in this situation, we can just safely remove the node from the tree, and we are done.
  B.        The node N to-be-removed has only one child. Here simply removing the node would disconnect the tree, but we can instead bypass it by connecting N’s parent to its children (independently of it being in the left or right subtree of N). This won’t violate any invariant for the BST: for instance, in the case shown in figure 9.15B, N is a left child of its parent P, so it’s smaller or equal than P, but likewise, all elements in the subtree rooted at N will be smaller or equal than P, including N’s child L.
  C.        If the node N that we need to remove has both children, we can’t just replace it with one of its children (for instance, if we were to replace the root in figure 9.15C with its right child R, then we would have to merge R’s left child with N’s, and that would require worst-case linear time (and also be a waste).

A.        The node we need to remove is a leaf: in this situation, we can just safely remove the node from the tree, and we are done.
  B.        The node N to-be-removed has only one child. Here simply removing the node would disconnect the tree, but we can instead bypass it by connecting N’s parent to its children (independently of it being in the left or right subtree of N). This won’t violate any invariant for the BST: for instance, in the case shown in figure 9.15B, N is a left child of its parent P, so it’s smaller or equal than P, but likewise, all elements in the subtree rooted at N will be smaller or equal than P, including N’s child L.
  C.        If the node N that we need to remove has both children, we can’t just replace it with one of its children (for instance, if we were to replace the root in figure 9.15C with its right child R, then we would have to merge R’s left child with N’s, and that would require worst-case linear time (and also be a waste). Instead, we can find the successor[46] of N: by construction, it will be the minimum node in its left subtree, let’s call it M, which in turn means the leftmost node of its right subtree. Once we have found it, we can delete M and replace N’s value with M’s value. No invariant will be violated, because by definition M will be no smaller than any node in N’s left branch (because it was no smaller than N, that was no smaller than any node in its left subtree), and no bigger than any node in N’s right branch (because it was its minimum).

(A)         Deleting a leaf.
(B)         Deleting a node with a single child. (Works symmetrically if the child is in the right branch).
(C)         Deleting a node with both children.

Figure 9.16 An example of how the successor of a node, or more in general the minimum of a subtree with respect to a certain coordinate, can lie anywhere in the subtree. The results of a few calls to findMin
on the whole tree and on node B
are explicitly shown.

For k-d trees, only leaves can be deleted easily: unfortunately, even if the node N
to be deleted has only one child C
, we can’t just replace N
with its child, because this would change the direction of the splits for C
and all of its subtree, as shown in figure 9.17.

function findMin(node, coordinateIndex)                                 #1
  if node == null then                                                  #2
    return null
  elsif node.level == coordinateIndex then                                      #3
    if node.left == null then                                           #4
      return node
    else                                                                        #5
      return findMin(node.left, coordinateIndex)
  else
    leftMin ← findMinNode(node.left, coordinateIndex)                  #6
    rightMin ← findMinNode(node.right, coordinateIndex)                #7
    return min(node, leftMin, rightMin)                                 #8
 

(A)        The initial tree, from which we remove node B. This node has only one child, so in a BST it would just be replaced by its child.
 (B)        However, in a k-d tree this might cause violating the k-d tree invariants, since moving a node one level up changes the coordinate on which the split is performed, when using it as a pivot.

If we replace N
the max of its left branch, however, it is possible that in N
’s old left branch there was another node with the same y
coordinate: in our example, that would be the case since there are two nodes with the same, maximal value for y
, node I
and node H
.
By moving node I
to replace B
, we would therefore break search, because node H
would never ever be found by the search
method in listing 9.3.

Luckily the solution is not too complicated: we can instead run findMin
on the left branch, replace N
’s point with the node M
found by findMin
, and set N
’s old left branch as the right branch of this new node we are creating, like we show in figure 9.19.
Then we just need to remove the old node M
from that right branch: notice that, differently from what happened with binary search trees, we can make no assumptions on M
here, so we might need to repeat these steps in the recursive call deleting M
.

Figure 9.19 Correct steps for deleting a node N
with only a left child. In figure 9.18 we have seen that finding the max of the left branch won’t work. Instead, we need to find the min M
, use it to replace the deleted node, and then set N
’s old left branch as the new node’s right branch. Then we only need to remove M
from the left branch, which requires a new call to remove (as you can see, unfortunately we can make no assumptions on this call, it might cascade and require more recursive calls to remove).
The worst case here is the latter: we don’t know where the minimum node will be in its branch, it could be far down the tree, or just the root of the branch; thus, if we end up in a  case similar to figure 9.17 (either with missing left or right branch, indifferently), we might have to call remove
recursively on n-1
nodes, in the absolute worst case.

function remove(node, point)                                            #1
  if node == null then                                                  #2
    return null
  elsif node.point == point then                                                #3
    if node.right != null then                                          #4
      minNode ← findMin(node.right, node.level)                        #5
      newRight ← remove(node.right, minNode.point)                     #6
      return new KdNode(minNode.point, node.left, newRight, node.level) #7
    elsif node.left != null then                                                #8
      minNode ← findMin(node.left, node.level)                         #9
      newRight ← remove(node.left, minNode.point)                      #10
      return new KdNode(minNode.point, null, newRight, node.level)      #11
    else
      return null                                                               #12
  elsif compare(point, node) < 0 then                                   #13
    node.left ← remove(node.left, point)
    return node
  else                                                          #14
    node.right ← remove(node.right, point)
    return node

Figure 9.20 The first few steps of nearest neighbor search. The first phase of nearest neighbor search consists of a search on the tree, during which we keep track of the distance of each node (aka point) we traverse: more precisely, if we need to find N
nearest neighbors for P
, we need to keep track of the N
smallest distances we find.
In this case we are showing the query for N=1
. Thus, if the search is successful, then we have for sure the nearest neighbor, at distance 0. Otherwise, when search ends, we are not yet sure we have found the actual nearest neighbor: in the example above, the point at the minimum distance we have found during traversal is D
but, as you can see from the figure, we can’t be sure that another branch of the tree doesn’t have a point within the radius given by dist(D)
.

How does nearest neighbor search work? We start from the consideration that each tree node covers one of the rectangular regions in which the space was partitioned, as we have shown in figures of 9.5 and 9.6. So first, we want to find which region contains our target point P
: that’s a good starting point for our search, because likely the point stored in the leaf covering that region, G
in this example, will be among the closest points to P
.
Can we be sure that G
will be the closest point to P
, though? That would have been ideal, but unfortunately that’s not the case. Figure 9.20 shows this first step in the algorithm, traversing a path in the tree from the root to a leaf, to find the smallest region containing P
.

function nearestNeighbor(node, target, (nnDist, nn)=(inf, null))        #1
  if node == null then                                                  #2
    return (nnDist, nn)
  else                                                          #3
    dist ← distance(node.point, target)                                #4
    if dist < nnDist then                                               #5
      (nnDist, nn) ← (dist, node.point)
    if compare(target, node) < 0 then                                   #6
      closeBranch ← node.left
      farBranch ← node.right
    else                                                                #7
      closeBranch ← node.right
      farBranch ← node.left
    (nnDist, nn) ← nearestNeighbor(closeBranch, target, (nnDist, nn)) #8
    if splitDistance(target, node) < nnDist then                        #9
      (nnDist, nn) ← nearestNeighbor(farBranch, target, (nnDist, nn))  #10
    return (nnDist, nn)                                                 #11

Figure 9.21 The second step, after an unsuccessful search, is to backtrack to check other branches of the tree for closer points. It is possible to have such points because when we traverse the tree, we cycle through the coordinates that we compare at each level, so we don’t always go in the direction of the closest point, but we are forced to go on one side of the pivot, depending only on a single coordinate: so it is possible that the nearest neighbor is on the wrong side of the pivot with respect to our target point P
.
In this example, when we reached D
, since it creates a vertical split, we needed to move to the left, as shown in figure 9.20: unfortunately the target point P
and its nearest neighbor lie on opposite sites of the split line for D
.

But even that is not enough: after finding the region containing P
, we can’t be sure that in neighboring regions there isn’t one or more points even closer than the closest point we have found during this first traversal.
Figure 9.21 exemplifies this situation perfectly: so far, we have found that D
is the closest point (among those visited) to P
, so the real nearest neighbor can’t be at a distance larger than the one between D
and P
: we can thus trace a circle (a hyper-sphere, in higher dimensions) centered at P
, and with radius equal to dist(D, P)
. If this circle intersects other partitions, then those regions might contain a point closer than D
, and we need to get to them.
How do we know if a region intersects our current nearest neighbor’s hyper-sphere? That’s simple: each region stems by the partitioning created by split lines; when traversing a path, we go on one side of the split (the one on the same side than P
), but if the distance between a split line and P
is less than the distance to our current nearest neighbor, then the hyper-sphere intersects the other partition as well.
To make sure we visit all partitions still intersecting the NN hyper-sphere, we need to back-track our traversal of the tree; in our example, we go back to node D
, then check the distance between P
and the vertical split line passing through D
(which, in turn, is just the difference of the x
coordinates of the two points), and since this is smaller than the distance to D
, our current NN, then we need to visit the other branch of D
as well. When we say visit here, we mean traversing the tree in a path from D
to the closest leaf to P
. While we do so, we visit node F
and we discover that it’s closer than D
, so we update our current NN (and its distance: you can see that we shrink the radius of our nearest neighbor perimeter, the circle marking the region where possible nearest neighbors can be found).
Are we done now? Not yet, we need to keep back-tracking to the root. We go back to node C
, but its split line is further away than our NN perimeter (and it didn’t have a right branch anyway), so we go back to node A
, as shown in figure 9.22.

Figure 9.22 We need to backtrack back towards the root. If we had to check every possible branch of the tree, however, this would be no better than scanning the whole list of points. Instead, we can use all the info we have to prune the search: in the geometric representation on the right, we show a visual hint of why it would be useless to check A
’s right sub-branch. If you check back figure 9.21, you can notice that, instead, we knew that we couldn’t rule out D
’s right sub-branch without traversing it.
At node A
we took the left branch during search, meaning P
is on the left semi-plane. If we were to traverse the right subtree of A
, all the points in that subtree would have their x
coordinate greater than or equal to A
’s. Therefore, the minimum distance between P
and any point in the right sub-tree of A
is at least the distance between P
and its projection on the vertical line passing through A
(i.e. A
’s split line): in other words, the minimum distance for any point right of A
will at least be the absolute value of the difference between P
’s and A
’s x
coordinates.
We can therefore prune search on A
’s right branch, and since it was the root, we are finally done. Notice how the more we climb back on the tree, the largest are the branches (and the regions) that we can prune – and so the largest is the saving.

·    Instead of the distance of a single point, we need to keep track of the m
shortest distances, if we want the m
closest points;
·    At each step, we use the distance of the m
-th closest point to draw our NN perimeter, and prune search;
·    To keep track of these m
distances, we can use a bounded priority queue. We have described something similar in chapter 2, section 2.7.3, when we described a method to find the m
largest values in a stream of numbers.
Listing 9.10 details the pseudo-code for the nNearestNeighbor
method.

function nNearestNeighbor(node, target, n)                              #1
  pq ← new BoundedPriorityQueue(n)                                     #2
  pq.insert((inf, null))                                                #3
  pq ← nNearestNeighbor(node, target, pq)                              #4
  (nnnDist, _) ← pq.peek()                                             #5
  if nnnDist == inf then                                                #6
    pq.top()
  return pq                                                             #7
 
function nNearestNeighbor(node, target, pq)                             #8
  if node == null then                                                  #9
    return pq
  else
    dist ← distance(node.point, target)                                #10
    pq.insert((dist, node.point))                                       #11
    if compare(target, node) < 0 then                                   #12
      closeBranch ← node.left
      farBranch ← node.right
    else                                                                #13
      closeBranch ← node.right
      farBranch ← node.left
    pq ← nNearestNeighbor(closeBranch, target, pq)                     #14
    (nnnDist, _) ← pq.peek()                                           #15
    if splitDistance(target, node) < nnnDist then                       #16
      pq ← nNearestNeighbor(farBranch, target, pq)                     #17
    return pq                                                           #18
 

(A)       An example built on the k-d tree from figures 9.20-9.22: by carefully choosing the target point, we can force the algorithm to search the whole tree – as shown on the tree representation.
  (B)       We show the spatial representation only of an edge case where all the points lie in a circle, and we choose the center of the circle as target of the search. Notice how the distance from P
to a split line will always be shorter than the radius of the circle (which, in turn, is the distance to the nearest neighbor).
So, there isn’t unfortunately much to do: the worst-case running time for this kind of query is O(n)
.

Nevertheless, to give you an intuition of why it works this way, consider a two-dimensional space: one point divides it into two halves, two points in three regions, 3 points will create four regions etc…, in general n
points will create n+1
regions; if the tree is balanced and n
is sufficiently big, we can assume these regions are approximately equally-sized.

In figure 9.25 you can see how we are going to prune branches: when we are at a split, we will certainly have to traverse the branch on the same side of the center of the search region P
; for the other branch, we check the distance between P
and its projection on the split line: if that distance is lower than or equal to the search region’s radius, it means that there is still an area of intersection between that branch and the search region, and so we need to traverse the branch; otherwise, we can prune it.

function pointsInSphere(node, center, radius)                           #1
  if node == null then                                                  #2
    return []
  else
    points ← []                                                        #3
    dist ← distance(node.point, center)                                #4
    if dist < radius then                                               #5
      points.insert(node.point)
    if compare(target, node) < 0 then                                   #6
      closeBranch ← node.left
      farBranch ← node.right
    else                                                                #7
      closeBranch ← node.right
      farBranch ← node.left
    points.insertAll(pointsInSphere(closeBranch, center, radius))       #8
    if splitDistance(target, node) < radius then                        #9
      points.insertAll(pointsInSphere(farBranch, center, radius))       #10
    return points                                                       #11

And similarly for vertical splits, checking the x
coordinates: likewise, it can be generalized for k-dimensional spaces, cycling through dimensions.

Figure 9.26 Region Search on a k-d tree: returns all the points in a k-d tree within a given hyper-rectangle. This means looking for points that, for each coordinate, satisfy two inequalities: each coordinate must be within a range; for instance, in this 2-D example, points’ x
coordinates need to be between -2
and 3.5
, and y
coordinates between -4
and 2
.

·    For very small regions intersecting only a single branch corresponding to a leaf, we will prune every other branch, and just follow a path to the leaf: so the running time will be O(h)
, where h
is the height of the tree – O(log(n))
for a balanced tree with n
points.
·    When the region is large enough to intersect all points, then we will have to traverse the whole tree, and the running time will be O(n)
.

function pointsInRectangle(node, rectangle)                             #1
  if node == null then                                                  #2
    return []
  else                                                          #3
    points ← []                                                        #4
    if rectangle[i].min≤node.point[i]≤rectangle[i].max " 0≤i<k then     #5
      points.insert(node.point)
    if intersectLeft(rectangle, node) then                              #6
      points.insertAll(pointsInRectangle(node.left, rectangle))
    if intersectRight(rectangle, node) then                             #7
      points.insertAll(pointsInRectangle(node.right, rectangle))
    return points                                                               #8

In high-dimensional spaces, the exponential term on k
for nearest neighbor becomes dominant, and makes not-worthy supporting the extra complexity of such a data structure.
Table 9.1 Operations provided by k-d tree, and their cost on a balanced k-d tree with n
elements

(*) Amortized, for a k-d tree holding k-dimensional points.
(§) When searching for the m
nearest neighbors, with m
constant (not a function of n
).

·   K-d trees work better on static datasets, because we can build balanced trees on construction, but insert
and remove
are not self-balancing operations.