Question: I have about 50,000 data points in 3D on which I have run scipy.spatial.Delaunay from the new scipy (I'm using 0.10) which gives me a very useful triangulation.

Based on: <http://en.wikipedia.org/wiki/Delaunay_triangulation> (section "Relationship with the Voronoi diagram")

...I was wondering if there is an easy way to get to the "dual graph" of this triangulation, which is the Voronoi Tesselation.

Any clues? My searching around on this seems to show no pre-built in scipy functions, which I find almost strange!

Thanks,
Edward

Answer: 
As I spent a considerable amount of time on this, I'd like to share my solution on how to get the Voronoi *polygons* instead of just the edges.

The code is at <https://gist.github.com/letmaik/8803860> and extends on the solution of [tauran](https://stackoverflow.com/a/15783581/60982).

First, I changed the code to give me vertices and (pairs of) indices (=edges) separately, as many calculations can be simplified when working on indices instead of point coordinates.

Then, in the `voronoi_cell_lines` method I determine which edges belong to which cells. For that I use the proposed solution of [Alink](https://stackoverflow.com/a/9471601/60982) from a related question. That is, for each edge find the two nearest input points (=cells) and create a mapping from that.

The last step is to create the actual polygons (see `voronoi_polygons` method). First, the outer cells which have dangling edges need to be closed. This is as simple as looking through all edges and checking which ones have only one neighboring edge. There can be either zero or two such edges. In case of two, I then connect these by introducing an additional edge. 

Finally, the unordered edges in each cell need to be put into the right order to derive a polygon from them.

The usage is:

```
P = np.random.random((100,2))

fig = plt.figure(figsize=(4.5,4.5))
axes = plt.subplot(1,1,1)

plt.axis([-0.05,1.05,-0.05,1.05])

vertices, lineIndices = voronoi(P)        
cells = voronoi_cell_lines(P, vertices, lineIndices)
polys = voronoi_polygons(cells)

for pIdx, polyIndices in polys.items():
    poly = vertices[np.asarray(polyIndices)]
    p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
    axes.add_patch(p)

X,Y = P[:,0],P[:,1]
plt.scatter(X, Y, marker='.', zorder=2)

plt.axis([-0.05,1.05,-0.05,1.05])
plt.show()

```

which outputs:

![Voronoi polygons](https://i.stack.imgur.com/LqND9.png)

The code is probably not suitable for large numbers of input points and can be improved in some areas. Nevertheless, it may be helpful to others who have similar problems.
I do not know of a function to do this, but it does not seem like an overly complicated task.

The Voronoi graph is the junction of the circumcircles, as described in the wikipedia article.

So you could start with a function that finds the center of the circumcircles of a triangle, which is basic mathematics (http://en.wikipedia.org/wiki/Circumscribed\_circle).

Then, just join centers of adjacent triangles.