Finding near neighbors

Question

I need to find "near" neighbors among a set of points.

There are 10 points in the above image. Red lines are edges from the Delaunay Triangulation, black stars mark the mid-lines of the edges, blue lines are the Voronoi tesselation. Point 1 has three "near" neighbors, i.e. 4, 6, and 7, but not 2 and 3, who are almost in line with the edge 1-7, but much further away.

What is a good way to identify the near neighbors (or "good" edges)? Looking at the figure, it seems to me that either selecting edges whose mid-point falls onto the intersection with the Voronoi lines, or considering as "near" neighbors those with touching Voronoi cells could be a good solution (the classification of 3-5 can go either way). Is there an efficient way of implementing either of the solutions in Matlab (I'd be happy to get a good general algorithm that I can then translate to Matlab, btw)?

Answer 1

You can implement your first idea of selecting edges whose mid-points fall on the intersection with the Voronoi lines by making use of the DelaunayTri class and its edges and nearestNeighbor methods. Here's an example with 10 random pairs of x and y values:

x = rand(10,1);                     %# Random x data
y = rand(10,1);                     %# Random y data
dt = DelaunayTri(x,y);              %# Compute the Delaunay triangulation
edgeIndex = edges(dt);              %# Triangulation edge indices
midpts = [mean(x(edgeIndex),2) ...  %# Triangulation edge midpoints
          mean(y(edgeIndex),2)];
nearIndex = nearestNeighbor(dt,midpts);  %# Find the vertex nearest the midpoints
keepIndex = (nearIndex == edgeIndex(:,1)) | ...  %# Find the edges where the
            (nearIndex == edgeIndex(:,2));       %#   midpoint is not closer to
                                                 %#   another vertex than it is
                                                 %#   to one of its end vertices
edgeIndex = edgeIndex(keepIndex,:);      %# The "good" edges

And now edgeIndex is an N-by-2 matrix where each row contains the indices into x and y for one edge that defines a "near" connection. The following plot illustrates the Delaunay triangulation (red lines), Voronoi diagram (blue lines), midpoints of the triangulation edges (black asterisks), and the "good" edges that remain in edgeIndex (thick red lines):

triplot(dt,'r');  %# Plot the Delaunay triangulation
hold on;          %# Add to the plot
plot(x(edgeIndex).',y(edgeIndex).','r-','LineWidth',3);  %# Plot the "good" edges
voronoi(dt,'b');  %# Plot the Voronoi diagram
plot(midpts(:,1),midpts(:,2),'k*');  %# Plot the triangulation edge midpoints

How it works...

The Voronoi diagram is comprised of a series of Voronoi polygons, or cells. In the above image, each cell represents the region around a given triangulation vertex which encloses all the points in space that are closer to that vertex than any other vertex. As a result of this, when you have 2 vertices that aren't close to any other vertices (like vertices 6 and 8 in your image) then the midpoint of the line joining those vertices falls on the separating line between the Voronoi cells for the vertices.

However, when there is a third vertex that is close to the line joining 2 given vertices then the Voronoi cell for the third vertex may extend between the 2 given vertices, crossing the line joining them and enclosing that lines midpoint. This third vertex can therefore be considered a "nearer" neighbor to the 2 given vertices than the 2 vertices are to each other. In your image, the Voronoi cell for vertex 7 extends into the region between vertices 1 and 2 (and 1 and 3), so vertex 7 is considered a nearer neighbor to vertex 1 than vertex 2 (or 3) is.

In some cases, this algorithm may not consider two vertices as "near" neighbors even though their Voronoi cells touch. Vertices 3 and 5 in your image are an example of this, where vertex 2 is considered a nearer neighbor to vertices 3 or 5 than vertices 3 or 5 are to each other.

Answer 2

I would agree that shared cell edges is a good neighbor criterion (based on this example). If you were using a mesh-oriented data structure (like something from Gts), then identifying shared edges would be trivial.

Matlab, on the other hand, makes this more "interesting". Assuming the voronoi cells are stored as patches, you might try obtaining the 'Faces' patch property (see this reference). That should return something like an adjacency matrix that identifies connected vertices. From that (and a little magic), you should be able to determine shared vertices, and then infer shared edges. In my experience, this sort of "search" problem is not well suited to Matlab - if possible, I recommend moving to a system more suited to queries of mesh connectivity.

Answer 3

Another possibility, simpler than creating triangulations or voronoi diagrams, is using a neighborhood matrix.

First place all your points into a 2D square matrix. Then you can run a full or partial spatial sort, so points will became ordered inside the matrix.

Points with small Y could move to the top rows of the matrix, and likewise, points with large Y would go to the bottom rows. The same will happen with points with small X coordinates, that should move to the columns on the left. And symmetrically, points with large X value will go to the right columns.

After you did the spatial sort (there are many ways to achieve this, both by serial or parallel algorithms) you can lookup the nearest points of a given point P by just visiting the adjacent cells where point P is actually stored in the neighborhood matrix.

You can read more details for this idea in the following paper (you will find PDF copies of it online): Supermassive Crowd Simulation on GPU based on Emergent Behavior.

The sorting step gives you interesting choices. You can use just the even-odd transposition sort described in the paper, which is very simple to implement (even in CUDA). If you run just one pass of this, it will give you a partial sort, which can be already useful if your matrix is near-sorted. That is, if your points move slowly, it will save you a lot of computation.

If you need a full sort, you can run such even-odd transposition pass several times (as described in the following Wikipedia page):

http://en.wikipedia.org/wiki/Odd%E2%80%93even_sort