# Best algorithm for calculating lists of neighbours

Given a collection of thousands of points in 3D, I need to get the list of neighbours for each particle that fall inside some cutoff value (in terms of euclidean distance), and if possible, sorted from nearest fo farthest.

Which is the fastest algorithm for this purpose?

PS: question also crossposted on TCS and SO

-
en.wikipedia.org/wiki/K-d_tree – whuber Jul 18 '11 at 17:25
corssposted on cstheory – Kaveh Jul 19 '11 at 4:19

Brute force is to bin (order n) then efficiently sort the bins with an n log n sort. This can even be done in parallel with parallel primitives (and runs fast on a GPU too). At that point you have to sort particles within bins and taking into account neighboring bins.

This is not optimal if you need to resort every timestep, particles only move one bin per timestep (most particle based codes) and or if the there is large range density range (in which case you may be sorting a large number of particles within a single bin).

When particles only move a small amount and you have a presorted system then an optimal resort is the way to go and I think this would be order N per spatial dimension.

If there is a wide range of density in your particle distribution you might consider a octal tree algorithm which is also N log N to make but would give you a nearest neighbor list based on the tree itself. There are efficient ways of fixing or updating trees as a particles system evolves. Trees can be constructed in parallel but I am not convinced they are very efficient yet, I think they use a Morton/hash to bin that also gives an octal tree structure and then they define tree nodes based on the bit changes in the bin ids. I have only seen implementations that go a few nodes deep in the tree, so there are ways to improve the parallel algorithms for this setting.

If you use the Morton hash to make your bins then the resulting sort can be used to count numbers of particles in bins, and bins twice or four times or 8 times as big. This is kind of like using the octal tree. Binning and sorting is pretty fast so you can do it a bunch of times depending upon what range of density you expect for your particle distribution.

Other ideas include using space filling curves to hash and other types of trees. The KD tree is popular but the parallel implementations of it I have seen have not been particularly clean or pretty.

-
what do you mean exactly by space fillling curves? can you point to some example? – flow Jul 19 '11 at 1:14
@flow the idea is to map a 3d or 2d space to a one d space using a curve for example en.wikipedia.org/wiki/Z-order_curve – Alice Jul 19 '11 at 1:36
this is hard to see/understand for me; some paper with a concrete application for this? – flow Jul 19 '11 at 3:04
If you are interested in GPU implementations I would read articles from the GPU GEMS series of books. The recent books have trees in them. GPU GEMS 3 had a nice description of a billiard simulation. The SDK implementation used a Morton hash that is also used by some of the trees. The space filling curves have the property that locality is preserved so the distance between two points on the 1 D curve is related to the distance in 3D. – Alice Jul 19 '11 at 12:09

You're looking for spherical range searching algorithms. I've had some success with k-D trees in lower dimensions.

-
which specific form of the algorithms did you use? for which application problem? – flow Jul 19 '11 at 0:58
@flow: Range searching on a k-D tree built from SIFT descriptors. – Jacob Jul 19 '11 at 14:29

Not sure if this is the best way, but to give you an idea on what you could do:

Sort on x, all points < distance into formula again (nlogn) Sort on y, all points < distance into formula (potential nlogn if you're unlucky) Compare all points < distance away and add them pairwise (if a is a neighbor to b then the reverse is true as well, on that note you should be able to reduce the amount of comparisons needed by removing points from the list which you are interested in..) Sort on distance once neighbours are complete... (nlogn again if you're unlucky).

-
sorry I do not understand your first sentence "sort on x, all points < distance into formula again" – flow Jul 19 '11 at 1:27
Use a standard sorting algorithm (nlogn) to sort all points by their x-value. Add all points within the cutoff in x-distance to a new vector and sort all of those on y. Now add those that are within the cutoff in y-distance. Sort the remaining by z and compare them. Add pairwise and make sure to flag it to not do double work. (if p0 is within distance of p1, then the reverse is true). It can be costly if you have a large density but could give you somewhere to start... – Ragshi Jul 19 '11 at 19:54

Not a general algorithm, but one recipe that can help (specially if your points move with time and you need to recompute your neighbour lists) is to divide the available space in cubes ("cells") of size $L$ (cutoff distance) (or some slightly larger value). Then you precompute/maintain, for each point, the cube to which it belongs, with a (perhaps bidirectional) mapping of "cube coordinates" => list of points inside . Then, when computing the neighbours of a point, you need only to consider the cube for this point and its eight neigbour cubes.

This is used in Molecular Dynamics simulations ("linked cells"), along with other related optimizations that you might find useful.

-
the problem i see with this (please correct me if I am wrong) is that dividing the space into the cubes and calculating to each point to which cube it belongs, it can take a lot of computation, or there is any fast way to do that? – flow Jul 19 '11 at 1:30
and furthermore, I only need to calculate the neighbours once, so there are not several steps like in molecular dynamics – flow Jul 19 '11 at 1:31
it all depends on the numbers involved. but to compute the cube for each point, is like (actually simpler) that computing one point-point distance, hence the total computation is of order $n$. OTOH, to compute all the raw pair distances is of order $n^2$ – leonbloy Jul 19 '11 at 1:44
yes, I got it now, thanks – flow Jul 19 '11 at 3:03

If you have enough points it may be worthwhile to compute the Delaunay triangulation of the points and use that to compute the neighbors. From the Delaunay triangulation finding the neighbors within a certain distance is a depth first search (include the neighbors within the bound, their neighbors within the bound etc.) In 3D the Delauany triangulation is difficult to program, but still fast, $O(n lg n)$, but there are libraries like CGAL that do it for you.

-
does this cgal library exist for gpu? – flow Jul 19 '11 at 1:14
@flow I do not know. I suspect that there are enough branch points that there is not enouh parallelism left for efficient gpu implementation, but I know very little about gpu coding. You might look at this forums.nvidia.com/index.php?showtopic=107315 which seems to be asking the same question. – deinst Jul 19 '11 at 1:28