4
$\begingroup$

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

  • 0
    [corssposted on cstheory](http://math.stackexchange.com/questions/52202/best-algorithm-for-calculating-lists-of-neighbours)2011-07-19

5 Answers 5

4

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.

  • 0
    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.2011-07-19
3

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

  • 0
    @flow: Range searching on a k-D tree built from SIFT descriptors.2011-07-19
2

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).

  • 0
    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...2011-07-19
2

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.

  • 0
    yes, I got it now, thanks2011-07-19
1

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.

  • 1
    @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 http://forums.nvidia.com/index.php?showtopic=107315 which seems to be asking the same question.2011-07-19