Nearest neighbor on a unit sphere, with roughly evenly distributed points.

Question

I'm writing a program that implements SCVT (Spherical Centroidal Voronoi Tesselation). I start with a set of points distributed over the unit sphere (I have an option for random points or an equal-area spiral). There will be from a several hundred to maybe 64K points.

I then need to produce probably several million random sample points, for each sample find the nearest point in the set, and use that to calculate a "weight" for that point. (This weigh may have to be looked up from another spherical set, but that set will stay static for any given run of the algorithm.)

Then I move the original points to the calculated points, and iterate the process, probably 10 or 20 times. This will give me the centers of the Voronoi tiles for subsequent use.

Later I will need to find a given point's nearest neighbor, to see what tile the user clicked on. This is trivially solved within the above problem, and doesn't need to be super-fast anyway. The part I need to be efficient is all those millions of nearest neighbors on the unit sphere. Any pointers?

Oh, I'm using x, y, z coordinates, but that's not set in stone. It just looks like it will simplify things. I'm also using C as I'm most familiar with it, but not wedded to that choice either. :)

I've considered using the spiral pattern for the sample points, as that gives me at least the last point's found neighbor as a good starting point for the next search. But if I do that, it looks like it would make any sort of tree search useless.

edit: [I'm sorry, I thought I was clear with the title and tags. I can generate random points easily. The issue is the nearest neighbor search. What's an efficient algorithm when all the points are on the unit sphere?]

Answer 1

Here is the article on neighbor search: http://en.wikipedia.org/wiki/Nearest_neighbor_search In my understanding you can use trivial algorithm of going through all Voronoi centers and calculate 3d distance between your point and center point.

distance_2 = (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2

where (x_0, y_0, z_0) is the point of interest (click) for you and {(x, y, z)} are Voronoi centers. The smallest distance will give you the nearest center.

Answer 2

Using a KD Trie is a good way to speed up the search. You can also get significantly better performance if you can tolerate some error. The ANN library will give you the result within an ε of your choosing.

Answer 3

OK. NEARPT3 http://www.ecse.rpi.edu/Homepages/wrf/Research/nearpt3/nearpt3.pdf algorithm could be helpful in your case. And it all depends on how many space you can afford to use for your N points. If it is O(N*logN) then there are algorithms like kD-tree based (http://www.inf.ed.ac.uk/teaching/courses/inf2b/learnnotes/inf2b-learn06-lec.pdf) which would work for O(logN) to find nearest point. In case of 64K point Nlog_2 N = about 10^6 which is easily can fit into memory of modern computer.

Answer 4

You may find that organising your points into a data structure called an Octree is useful for efficient search for nearby points. See http://en.wikipedia.org/wiki/Octree

Answer 5

Your points are uniformly distributed over the sphere. Therefore, it would make a lot of sense to convert them to spherical coordinates and discretize. Searching the 2D grid first would narrow down the choice of nearest neighbour to a small part of the sphere in constant time.

Answer 6

I have devised a curve (I'm sure I'm not the 1st) that spirals along the sphere from pole to pole. It remains a constant distance from neighboring windings (if I did it right). For z [-1 at south pole to +1 at north pole]:

n= a constant defining a given spiral k=sqrt(n*pi)

r=sqrt(z^2) theta=k*asin(z) x=r*cos(theta) y=r*sin(theta)

It makes k/2 revolutions around the sphere, with each winding sqrt(4pi/n) from adjacent windings, while the slope dz/d(x,y) is 1/k.

Anyway, set k such that the inter-winding distance covers the largest tile on the sphere. For every point in the main set, calculate the theta of the nearest point on the curve, and index the list of points by those numbers. For a given test point, calculate it's (theta of the nearest point on the curve), and find that in the index. Search outward (in both directions) from there, to theta values that are as far away as your current nearest neighbor. After reaching that limit, if the distance to that neighbor is less than the distance from the test point to the next adjacent winding, you've found the nearest neighbor. If not, jump the theta value by 2pi and search that winding the same way.

Critique?