KD Tree - Nearest Neighbor Algorithm

Question

I'm not quite understanding the O(log n) nearest neighbor algorithm from wikipedia.

…

…

The algorithm unwinds the recursion of the tree, performing the following steps at each node:

...

The algorithm checks whether there could be any points on the other side of the splitting plane that are closer to the search point than the current best. In concept, this is done by intersecting the splitting hyperplane with a hypersphere around the search point that has a radius equal to the current nearest distance. Since the hyperplanes are all axis-aligned this is implemented as a simple comparison to see whether the difference between the splitting coordinate of the search point and current node is less than the distance (overall coordinates) from the search point to the current best.

If the hypersphere crosses the plane, there could be nearer points on the other side of the plane, so the algorithm must move down the other branch of the tree from the current node looking for closer points, following the same recursive process as the entire search.

If the hypersphere doesn't intersect the splitting plane, then the algorithm continues walking up the tree, and the entire branch on the other side of that node is eliminated.

It's 3.2 that is confusing me and I have seen this question. I'm implementing the algorithm in Java and not sure if I got it right.

//Search children branches, if axis aligned distance is less than current distance
if (node.lesser!=null) {
    KdNode lesser = node.lesser;
    int axis = lesser.depth % lesser.k;
    double axisAlignedDistance = Double.MAX_VALUE;
    if (axis==X_AXIS) axisAlignedDistance = Math.abs(lastNode.id.x-lesser.id.x);
    if (axis==Y_AXIS) axisAlignedDistance = Math.abs(lastNode.id.y-lesser.id.y);
    else if (axis==Z_AXIS) axisAlignedDistance = Math.abs(lastNode.id.z-lesser.id.z);

    //Continue down lesser branch
    if (axisAlignedDistance<=lastDistance && !set.contains(lesser)) {
        searchNode(value,lesser,set,K);
    }
}
if (node.greater!=null) {
    KdNode greater = node.greater;
    int axis = greater.depth % greater.k;
    double axisAlignedDistance = Double.MAX_VALUE;
    if (axis==X_AXIS) axisAlignedDistance = Math.abs(lastNode.id.x-greater.id.x);
    if (axis==Y_AXIS) axisAlignedDistance = Math.abs(lastNode.id.y-greater.id.y);
    else if (axis==Z_AXIS)axisAlignedDistance = Math.abs(lastNode.id.z-greater.id.z);

    //Continue down greater branch
    if (axisAlignedDistance<=lastDistance && !set.contains(greater)) {
        searchNode(value,greater,set,K);
    }
}

Does the code above accomplishing the 3.2 aspect of the algorithm? Specifically where I populate the "axisAlignedDistance" variable.

You can find the full source code of the KDTree here.

Thanks for any help/pointers.

have you tested the code? Why ask a question about code that you think will do what it's supposed to? Without testing the code myself, i'd say it looks like it should do the job. And, i would avoid while(true) loops if possible. Especially because in your case, after the loop it is always current = node. — kutschkem, Jun 26 '12 at 14:23
Better to test the code and then ask questions about understanding the algorithm. IF it doesn't work, then you can post output and we can go from there. — Chris Dargis, Jun 26 '12 at 14:46
@kutschkem I've tested the code and it returns the correct results but that doesn't mean it implemented the algorithm correctly. I mean a O(n) algorithm would return the correct results but would not be correct. I'll update the question to help narrow down the question. — Justin, Jun 26 '12 at 16:28
@VanDarg I've tested the code, I'll also update the question to help narrow down the actual question. — Justin, Jun 26 '12 at 16:29
@Justin: Sounds good. It may also be a good idea to specify the Big-O bound you are looking for ( if O(N) is not sufficient ) — Chris Dargis, Jun 26 '12 at 16:31
@Justin: Are you confused if this algorithm runs in O(log(N)) time? Is the method written in your question the **complete** searchNode() function? Of what type is the `set` object? Is it possible to also see the `contains` method? — Chris Dargis, Jun 26 '12 at 16:39
@VanDarg I've linked the full source code in the question. I'm confused if my algorithm is following the algorithm in question. The algorithm should avoid branches of the tree according to 3.2 of the NNS algorithm. I am just not sure if my code is following 3.2 since I don't quite understand that part of it. — Justin, Jun 26 '12 at 16:45
I found a problem with the original code. I've updated my question to address the problems I found. — Justin, Jun 27 '12 at 17:36

Justin · Accepted Answer · 2013-11-21T14:00:56.580

I am adding this, hoping it helps others that have search out the same problem. I ended up addressing 3.2 with the following code. Although, I am not sure if it is 100% correct. It has passed all the tests I have come up with. The original code above has failed on a number of the same test cases.

More explicit solution using Point, Line, Rectangle, and Cube objects:

int axis = node.depth % node.k;
KdNode lesser = node.lesser;
KdNode greater = node.greater;

//Search children branches, if axis aligned distance is less than current distance
if (lesser!=null && !examined.contains(lesser)) {
    examined.add(lesser);

    boolean lineIntersectsRect = false;
    Line line = null;
    Cube cube = null;
    if (axis==X_AXIS) {
        line = new Line(new Point(value.x-lastDistance,value.y,value.z), new Point(value.x+lastDistance,value.y,value.z));
        Point tul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tur = new Point(node.id.x,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tlr = new Point(node.id.x,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Rectangle trect = new Rectangle(tul,tur,tlr,tll);
        Point bul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bur = new Point(node.id.x,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point blr = new Point(node.id.x,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Rectangle brect = new Rectangle(bul,bur,blr,bll);
        cube = new Cube(trect,brect);
        lineIntersectsRect = cube.inserects(line);
    } else if (axis==Y_AXIS) {
        line = new Line(new Point(value.x,value.y-lastDistance,value.z), new Point(value.x,value.y+lastDistance,value.z));
        Point tul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tlr = new Point(Double.POSITIVE_INFINITY,node.id.y,Double.NEGATIVE_INFINITY);
        Point tll = new Point(Double.NEGATIVE_INFINITY,node.id.y,Double.NEGATIVE_INFINITY);
        Rectangle trect = new Rectangle(tul,tur,tlr,tll);
        Point bul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point blr = new Point(Double.POSITIVE_INFINITY,node.id.y,Double.POSITIVE_INFINITY);
        Point bll = new Point(Double.NEGATIVE_INFINITY,node.id.y,Double.POSITIVE_INFINITY);
        Rectangle brect = new Rectangle(bul,bur,blr,bll);
        cube = new Cube(trect,brect);
        lineIntersectsRect = cube.inserects(line);
    } else {
        line = new Line(new Point(value.x,value.y,value.z-lastDistance), new Point(value.x,value.y,value.z+lastDistance));
        Point tul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tlr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Rectangle trect = new Rectangle(tul,tur,tlr,tll);
        Point bul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,node.id.z);
        Point bur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,node.id.z);
        Point blr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,node.id.z);
        Point bll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,node.id.z);
        Rectangle brect = new Rectangle(bul,bur,blr,bll);
        cube = new Cube(trect,brect);
        lineIntersectsRect = cube.inserects(line);
    }

    //Continue down lesser branch
    if (lineIntersectsRect) {
        searchNode(value,lesser,K,results,examined);
    }
}
if (greater!=null && !examined.contains(greater)) {
    examined.add(greater);

    boolean lineIntersectsRect = false;
    Line line = null;
    Cube cube = null;
    if (axis==X_AXIS) {
        line = new Line(new Point(value.x-lastDistance,value.y,value.z), new Point(value.x+lastDistance,value.y,value.z));
        Point tul = new Point(node.id.x,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tlr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tll = new Point(node.id.x,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Rectangle trect = new Rectangle(tul,tur,tlr,tll);
        Point bul = new Point(node.id.x,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point blr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bll = new Point(node.id.x,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Rectangle brect = new Rectangle(bul,bur,blr,bll);
        cube = new Cube(trect,brect);
        lineIntersectsRect = cube.inserects(line);
    } else if (axis==Y_AXIS) {
        line = new Line(new Point(value.x,value.y-lastDistance,value.z), new Point(value.x,value.y+lastDistance,value.z));
        Point tul = new Point(Double.NEGATIVE_INFINITY,node.id.y,Double.NEGATIVE_INFINITY);
        Point tur = new Point(Double.POSITIVE_INFINITY,node.id.y,Double.NEGATIVE_INFINITY);
        Point tlr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Point tll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY);
        Rectangle trect = new Rectangle(tul,tur,tlr,tll);
        Point bul = new Point(Double.NEGATIVE_INFINITY,node.id.y,Double.POSITIVE_INFINITY);
        Point bur = new Point(Double.POSITIVE_INFINITY,node.id.y,Double.POSITIVE_INFINITY);
        Point blr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Rectangle brect = new Rectangle(bul,bur,blr,bll);
        cube = new Cube(trect,brect);
        lineIntersectsRect = cube.inserects(line);
    } else {
        line = new Line(new Point(value.x,value.y,value.z-lastDistance), new Point(value.x,value.y,value.z+lastDistance));
        Point tul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,node.id.z);
        Point tur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,node.id.z);
        Point tlr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,node.id.z);
        Point tll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,node.id.z);
        Rectangle trect = new Rectangle(tul,tur,tlr,tll);
        Point bul = new Point(Double.NEGATIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bur = new Point(Double.POSITIVE_INFINITY,Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point blr = new Point(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Point bll = new Point(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY);
        Rectangle brect = new Rectangle(bul,bur,blr,bll);
        cube = new Cube(trect,brect);
        lineIntersectsRect = cube.inserects(line);
    }

    //Continue down greater branch
    if (lineIntersectsRect) {
        searchNode(value,greater,K,results,examined);
    }
}

I think this simpler code should also work, it has passed all the same tests as the above code.

int axis = node.depth % node.k;
KdNode lesser = node.lesser;
KdNode greater = node.greater;

//Search children branches, if axis aligned distance is less than current distance
if (lesser!=null && !examined.contains(lesser)) {
    examined.add(lesser);

    double p1 = Double.MIN_VALUE;
    double p2 = Double.MIN_VALUE;
    if (axis==X_AXIS) {
        p1 = node.id.x;
        p2 = value.x-lastDistance;
    } else if (axis==Y_AXIS) {
        p1 = node.id.y;
        p2 = value.y-lastDistance;
    } else {
        p1 = node.id.z;
        p2 = value.z-lastDistance;
    }
    boolean lineIntersectsCube = ((p2<=p1)?true:false);

    //Continue down lesser branch
    if (lineIntersectsCube) {
        searchNode(value,lesser,K,results,examined);
    }
}
if (greater!=null && !examined.contains(greater)) {
    examined.add(greater);

    double p1 = Double.MIN_VALUE;
    double p2 = Double.MIN_VALUE;
    if (axis==X_AXIS) {
        p1 = node.id.x;
        p2 = value.x+lastDistance;
    } else if (axis==Y_AXIS) {
        p1 = node.id.y;
        p2 = value.y+lastDistance;
    } else {
        p1 = node.id.z;
        p2 = value.z+lastDistance;
    }
    boolean lineIntersectsCube = ((p2>=p1)?true:false);

    //Continue down greater branch
    if (lineIntersectsCube) {
        searchNode(value,greater,K,results,examined);
    }
}

Chris Dargis · Answer 2 · 2012-06-26T19:04:12.677

It is important to note the following from 3:

"The algorithm unwinds the recursion of the tree, performing the following steps at each node: "

Your implementation appears to be only making recursive calls, and not anything upon the unwinding of the recursion.

Although it appears you are detecting intersection correctly, you are not performing the steps upon the backing out (unwinding) of recursion. After your recursive calls are made there are 0 steps taken (including 3.2).

3.2 states: "If the hypersphere doesn't intersect the splitting plane, then the algorithm continues walking up the tree, and the entire branch on the other side of that node is eliminated"

What does this mean? After you reach a base case (algorithm reaches a leaf node), the recursion begins to unwind. After each level of recursion is unwound, the algorithm checks to see if a sub-tree could possibly contain a closer neighbor. If it can, another recursive call is made on that sub-tree, if not the algorithm continues to unwind (walks up the tree).

You should start with this:

"Starting with the root node, the algorithm moves down the tree recursively, in the same way that it would if the search point were being inserted"

And then start thinking about what steps to take and when, during the unraveling of the recursion.

This is a challenging algorithm. If you have completed the insert() method of this data-structure, you can use that as a beginning framework for this algorithm.

EDIT:

//Used to not re-examine nodes
Set<KdNode> examined = new HashSet<KdNode>();

It may be easier, faster, and would use less memory to simply place a flag in each node that you can mark as "visited". Ideally, HashSets have a constant time lookup, but there really is no way to guarantee this. That way, you don't have to search through a set of nodes upon each visit. Before doing a search of the tree, you can initialize all of these values to false in O(N) time.

The unwinding should occur in the "nearestNeighbourSearch()" method call which calls the "searchNode()" initially on the leaf node. After the "searchNode()" finishes on the leaf node, it returns to the "nearestNeighbourSearch()" method which moves to the leaf's parent and calls "searchNode()" on it. — Justin, Jun 26 '12 at 18:37
@Justin: I see that now. But you have two separate methods traversing the tree here. You need to be performing these checks after each level of recursion unwinds in order to ensure the run-time of `O(log(N))`. — Chris Dargis, Jun 26 '12 at 18:51
about the examined Set. I plan to do just that when I get the algorithm working correctly. — Justin, Jun 27 '12 at 19:49

score 0 · Answer 3 · edited May 23 '17 at 12:31

0

Yes, the description of NN (Nearest Neighbour) search in a KD Tree on Wikipedia is a little hard to follow. It doesn't help that a LOT of the top Google search results on NN KD Tree searches are just plain wrong!

See https://stackoverflow.com/a/37107030/591720 for a correct explanation/algorithm.

edited May 23 '17 at 12:31

Community

1
1

answered May 09 '16 at 02:39

Scott Smedley

1,779
20
28

KD Tree - Nearest Neighbor Algorithm

3 Answers3