QuadTree find neighbor

There are some known algorithms. Check them out.

1. Kunio Aizawa et al. - Constant Time Neighbor Finding in Quadtrees: An Experimental Result
2. Kasturi Varadarajan - All Nearest Neighbours via Quadtrees
3. Robert Yoder, Peter Bloniarz - A Practical Algorithm for Computing Neighbors in Quadtrees, Octrees, and Hyperoctrees

If your language has good support for arrays and you can choose the representation of your tree, then this turns out to be much simpler than various papers suggest.

Trick is to represent the parent-child relationship as a vector:

def neighbour(tree, node, direction):
    """Finds the neighbour of a node in an adaptive quadtree or it's D-dimensional
    generalization (orthantree?).

    Be very careful with indexing when implementing this. Because it indexes  
    physical space with arrays, sometimes you need to think of coords in terms
    of Euclidean coords, and other times in terms of array coords.

    Args:
        tree: an object holding a bunch of node attributes, indexed by node:
            * `parent`, an (N,)-array of integers. The `n`th element gives the 
               index of `n`'s parent. The parent of the root is -1.
            * `children`, an ((N,) + (2,)*D)-array of integers. The `n`th slice 
               gives the indices of `n`'s children. The bottom-left child of node 3
               in a quadtree would be (3, 0, 0). 
            * `descent`, an (N, D)-array with elements from {-1, +1}. The `n`th
               row gives which direction node `n` lies in compared to its parent.
               For example, the left-bottom quadrant in a quadtree would be `(-1, -1)`.
            * `terminal`, an (N,)-array of booleans. The `n`th element is True
               if node `n` is a leaf.
        node: an integer, the index of the node you want to find the neighbour of.
        direction: a (D,)-array with elements from {-1, +1}

    Returns:
        An integer giving the index of the neighbouring node, or -1 if it doesn't 
        exist.
    """
    direction = np.asarray(direction)

    # Ascend to the common ancestor
    neighbour_descents = []
    while True:
        if (direction == 0).all() or node < 0:
            break
        node_descent = tree.descent[node]
        neighbour_descent = node_descent*(1 - 2*abs(direction))
        neighbour_descents.append(neighbour_descent)

        direction = ((node_descent + direction)/2).astype(int)
        node = tree.parent[node]

    # Descend to the neighbour 
    for neighbour_descent in neighbour_descents[::-1]:
        if tree.terminal[node] or node < 0:
            break
        node = tree.children[(node, *(neighbour_descent.T + 1)//2)]

    return node

It supports bitrees (?), quadtrees, octrees, and general N-dimensional trees (hyperoctree? orthantree?). It also supports any direction - cardinal or diagonal. Finally, it's really easy to vectorize.

Inspiration was the FSM-based approach from Yoder that @torvin posted, plus a requirement this work for any number of dimensions.

There's test and demo code here.