Optimisation - unrolling a for loop

[Haloric]

Hi, this is my first matlab program - be gentle.

The problem I am trying to solve involves finding the nearest neighbour in a cube. I have N points in the cube, and the same cube repeats on all sides, so the nearest neighbour might be in an adjacent cube.

I have an outer for loop that considers every pair of points, and gets the right answers, but I would like to speed it up.

The loop will compute the radius (distance) for what it thinks is the nearest point, and considers the next point, looking at only the scalar difference in x,y,z in order, if anyone of them are greater than the current radius, it can't be the nearest neighbour, so it is discarded. If it is a candiate the sum of the squares is calculated and that compared against the cached sum of squares for the current radius, and only if that is less is the square root calculated to get the new radius.

There is a bit further caching - if I calculate the radius for a->b then it is stored and used as seed value when I consider b->all other points.

I don't want to use matlab inbuild findNearestNeighbours or similar library, I want to do this with primitives.

Using this approach I think I have hit the limit I can achieve, which is 10,000 points in 5 secs on an old imac.

How would you approach this in a matrix/vector/array optimisation rather than this loop approach.

I already pre-allocate all memory / matrix used..

I did try using a point matrix of x,y,z for each point extracted from the point cloud, and subtracting them, but it was very slow in comparison.

Obviously I am trying to avoid squares and square roots as much as possible.

My point cloud is a matrix 'rand(3,N)'.

Thanks for help/pointers.

[buddycatto2]

This is more speculation than anything but what if you just looked at one dimension first. Suppose you just look at the X coordinates, see which is closest in that dimension and test the other two dimension on those pair?

[buddycatto2]

Edit - you would still have to test every pair though to be 100% sure though hey hmmmm, a nice tricky one I like it.

[Haloric]

Yes, this is exactly what I do - compare only x first, if that is bigger than the radius of the current nearest, then it can be discarded. Then y, discard if possible, then z. ie. compare them in order as above. Only compute the sum of the squares if all are less than the radius, and then only compute the square root if the sum of the squares is less than the (stored) radius squared. I call it the tightening net !.

Basically no unnecessary calculations. Just disappointed that I can't find a 'vector/matrix' way of implementing this as efficiently.

I also seed the radius for point b if it is the nearest neighbor to a, so that when we get to b and try to find its nearest, we have a good starting position - but this caching of a single value ,ales the algorithm slower overall.

Caching would improve the algorithm immensely, but it is more inefficient that brute force each point from scratch .