How can I get a dictionary of cells from this Voronoi Diagram data?

Question

Using the voronoi/delaunay diagram generation library found in this program, which is based on Fortune's original implementation of his algorithm, with a random set of points as input data, I am able to get the following output data:

1. A list of the edges from the Delaunay Triangulation, meaning that for each input point, I can see which input points are its neighbors. They don't appear to be in any particular order.
2. A list of the vertex pairs from the Voronoi Diagram, which I can use to draw the Voronoi diagram one line at a time. Again, apparently in no particular order.
3. An unnamed list of pairs of points, which seems to just be the same list as 2, but in a different order.
4. A list of the vertices formed in the Voronoi Diagram, also apparently in no particular order.

Here is an example of data from a test run of my program using this library:

Input points:
0   (426.484, 175.16)
1   (282.004, 231.388)
2   (487.891, 353.996)
3   (50.8574, 5.02996)
4   (602.252, 288.418)

Vertex Pairs: 
0   (387.425, 288.533)  (277.142, 5.15565)
1   (387.425, 288.533)  (503.484, 248.682)
2   (277.142, 5.15565)  (0, 288.161)
3   (387.425, 288.533)  (272.213, 482)
4   (503.484, 248.682)  (637.275, 482)
5   (503.484, 248.682)  (642, 33.7153)
6   (277.142, 5.15565)  (279.477, 0)

Voronoi lines?: 
0   (279.477, 0)    (277.142, 5.15565)
1   (642, 33.7153)  (503.484, 248.682)
2   (503.484, 248.682)  (637.275, 482)
3   (387.425, 288.533)  (272.213, 482)
4   (277.142, 5.15565)  (0, 288.161)
5   (387.425, 288.533)  (503.484, 248.682)
6   (277.142, 5.15565)  (387.425, 288.533)

Delaunay Edges: 
0   (282.004, 231.388)  (487.891, 353.996)
1   (602.252, 288.418)  (487.891, 353.996)
2   (426.484, 175.16)   (487.891, 353.996)
3   (426.484, 175.16)   (602.252, 288.418)
4   (50.8574, 5.02996)  (282.004, 231.388)
5   (426.484, 175.16)   (282.004, 231.388)
6   (50.8574, 5.02996)  (426.484, 175.16)

Vertices: 
0   (277.142, 5.15565)
1   (503.484, 248.682)
2   (387.425, 288.533)
3   (0, 288.161)
4   (272.213, 482)
5   (637.275, 482)
6   (642, 33.7153)
7   (279.477, 0)

While the above data is adequate if all I need is to draw the Voronoi and Delaunay diagrams, it is not enough information for the actual work I am trying to do with these diagrams. What I need is a dictionary of polygons formed by the Voronoi vertices, indexed by the input point that each polygon was formed around. Preferably, for each polygon, these points would be sorted in clockwise order.

With the above information, I could implicitly assign data to each region, assign data to corners if necessary, tell which regions share edges (using the Delaunay edges), and do analysis accordingly.

So in short, how can I use the data available to me to put together a dictionary in which the key is one of the input points, and the data indexed by that key is a list of the Voronoi vertices that form the surrounding polygon? Or alternatively, is that information somewhere implicit in the data I've been given?

Insufflate

Is this all you get from the library? Some of the Voronoi cells are not described with a closed polygon. – Silurid

This is all the library gives me. The voronoi cells not described by a closed polygon are (I think) cells that meet the edge of the rectangular plane; for example, all of the border cells in this diagram: www-cs-students.stanford.edu/~amitp/game-programming/… – Insufflate

Answer 1

Fortune's algorithm is O(n log n) - but your code will be O(n^2), if you try to reconstruct cells brute-force fashion as proposed by Alink.

The starting point for my answer is that what you are using to generate the cells is not a library, but rather is just a class written to neatly wrap up the code originally presented by Fortune himself, and not actually a mature library. So, the author in fact hasn't anticipated your needs, and although the information you want has been computed, it isn't accessible.

Internally, your input points are stored as instances of the "Site" struct, and the algorithm proceeds to create half-edges, each of which maintains a reference "vertex" which is a pointer to the Site it encloses. Stepping along half-edges you naturally circumnavigate the enclosed Site - exactly what you need.

To access this data, I suggested modifying or extending the VoronoiDiagramGenerator class; I would do it by creating a hash table with Site pointers as the key and a single HalfEdge pointer as the value. Then, modify the generateVoroni method, inserting your new code immediately following the call to voronoi:

For each HalfEdge in ELHash
         Get table entry for current half edge's Site
         If site in table has null HalfEdge reference
            set current HalfEdge reference
         End If
End For each

...and there is your dictionary. That single half-edge will allow you to "walk" the perimeter of the polygon enclosing the related site, which I think is what you asked for. Your next problem will be to efficiently discover which polygon encloses some new data point - but that is another question :-). I hope you'll consider sharing your completed class - it should be a significantly more useful than the base class.

Edit: Here is an excellent presentation descibing all said above in pictures: http://ima.udg.es/~sellares/ComGeo/Vor2D_1.ppt:

• definition of Voronoy Diagram
• tree of of half-edges (see pics. below)
• Fortunes algorithm in pictures

And here is a C# implementation which could help you to retrieve the dictionary, as proposed above: http://www.codeproject.com/Articles/11275/Fortune-s-Voronoi-algorithm-implemented-in-C

Affectional

Well, I think my algorithm could be improved under O(N^2) by optimizing the nearest search with some space partitioning method (spiral in grid, quadtree ...). Such structure could also be useful if he needs to quickly find which voronoi cell contains a random point. But I agree with you, the library should provide this. – Haik

You are right there, I was thinking only of a very naive implementation. As you say, it would be be great for quickly finding enclosing cells. Be nice if this very useful code could be shared like the original class. – Affectional

Thanks for the great edit Denis - half-edges are definitely a lot clearer when illustrated. – Affectional

Just as I had done the relevant changes to the class, I found an already existing modification online: hep.wisc.edu/~dasu/public/Delphes-3.0.9/external/fastjet/… hep.wisc.edu/~dasu/public/Delphes-3.0.9/external/fastjet/… (With small changes, it is usable stand alone.) For each edge, it returns the indices of the two sites. – Paramagnetism

Answer 2

Your list of edges is somehow incomplete, you need to add the ones at the border of the containing rectangle provided to the library call (seems to be 642,482 here). Technically, a Voronoi subdivision should use some infinite edges, but those are all finite. I assume that you also want these "open" polygons near this border, since they are all like that in your example.

Adding those border edges seem not hard, just tedious. Probably something like, for each side of the main rectangle, find all vertices on it (ignoring corners), sort them (by x for the horizontal one, by y for vertical) and split that side using these values. This generates the missing edges, but don't add them directly to your main list, because they are special since they are the only ones not separating two cells.

So, for the question itself, I would go like this: In your main list of edges (provided by the library), each edge separates two cells and if we find which ones, then we can just assign that edge to each one of these cells. Since a cell is equivalent to an input point, we will have the dictionary wanted, except with a list of edges instead of vertices, but that's easy to convert.

Now to get these 2 cells: Calculate the middle point of the edge and from this, find the two nearest input points by simply iterating through the list while keeping the 2 smallest distances. By the properties of the Voronoi structure, those two are the ones forming the two cells. Note that these two distances should be equal, but the float imprecision will probably introduce a slight difference.

To finish, add the border edges that we generated along the main rectangle, but for those, just use the first nearest input point, since they are only adjacent to one cell.

Finally, we can convert each list of edges to a list of vertices (dump each ends points into a set). If you want to sort them in clockwise order, note that it's a convex polygon with an input point inside. So, you can just generate the vector going from the input point to each vertex, calculate its angle from one axis (use std::atan2(x,y)) and use this angle as comparator value to sort them (see std::sort).

Haik

Answer 3

I used Triangle package for generating Dalaunay triangulation: http://www.cs.cmu.edu/~quake/triangle.html

It works in 2 modes a) as a triangulate.exe utility and b) as a C library To compile it as a utility you just need to compile triangle.c and run:

   triangulate -vZ input.poly 
   #v -voronoy, Z - starting from 0 index

to get voronoi diagram (Refer to the manual about .poly format) I've made an experiment with your input data in such a .poly file:

# <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> 
5 2 0 0
# Following lines: <vertex #> <x> <y> [attributes] [boundary marker]
0 426.484 175.16
1 282.004 231.388
2 487.891 353.996
3 50.8574 5.02996
4 602.252 288.418
#One line: <# of segments> <# of boundary markers (0 or 1)>
5 0
#Following lines: <segment #> <endpoint> <endpoint> [boundary marker]
0 0 1
1 1 2
2 2 3
3 3 4
4 4 0

But it just reports input data error.

• Working with this package I'd say that it often doesn't work with with ïnput data reporting an error. It accepts only input polygons (not random points), and the issue here is that you have self-intersecting input polygon.
• It doesn't answer your question, reporting just set of points, not a dictionary

Lens

Answer 4

You could just use Triangle: http://www.cs.cmu.edu/~quake/triangle.html

Silurid

Questionable licensing and complete lack of a programming reference aside, after looking at triangle.h to try and see the available API, I don't really see how this solves my problem. Please clarify? – Insufflate

Answer 5

http://svn.osgeo.org/qgis/trunk/qgis/python/plugins/fTools/tools/voronoi.py

This python implementation by Carson Farmer gives you topological info

Scheck