Voronoi Diagram Construction

Joseph O’Rourke, 18 Jan 2012.
>>Use as you see fit, and give credit as seems appropriate.<<

This notebook illustrates two different methods for constructing the Voronoi Diagram in Mathematica.  The first uses the ComputationalGeometry package.  The second uses (or rather, misues!) the ListDensityPlot[ ] function.

We start with the ComputationalGeometry package:

Needs["ComputationalGeometry`"]

Generate n random points in a unit square [0,1]^2, for testing purposes.  Leave n global:

n=50;
pts=Table[{RandomReal[],RandomReal[]},{n}];

DiagramPlots[] with all default choices scales so that every Voronoi vertex is in the image, which is generally NOT desired.  For example:

DiagramPlot[pts]

Note: I am using a hack to color the diagrams.  I couldn’t figure out how to do it more directly.  So after the graphics is produced, I change every instance of a Line[] to a colored Line[] via a substitution.  (Anyone who knows a better way, please let me know!)

The option TrimPoints is used to remove distant Voronoi vertices that distort the scale.

g1=DiagramPlot[pts,LabelPoints->False,TrimPoints->5]/.Line[x_]->{Blue,Line[x]}

Now PlanarGraphPlot[ ] can be used to draw the Delaunay triangulation.  Again artificial coloring hack:

g2=PlanarGraphPlot[pts,LabelPoints->False]/.Line[x_]->{Red,Line[x]}

Display the two on top of one another:

Show[g1,g2]

Now I turn to the second method, using the new  ListDensityPlot[ ] function.

The options to RandomReal[] cause it to generate triples, the third number of which is used to control the colors via a ColorFunction.  The key reason this ListDensityPlot[ ] works is that InterpolationOrder is set to zero, which means that it fits constants (polynomials of degree 0) to interpolate.  Then the ListDensityPlot[ ] function decides to partition the regions according to the Voronoi diagram!  So the Voronoi diagram computation is hidden underneath the ListDensityPlot[ ] function.

pts3=RandomReal[1,{50,3}];g1=ListDensityPlot[pts3,InterpolationOrder->0,ColorFunction->"SouthwestColors",Mesh->All]

Now strip off the 3rd coordinate so the sites may be plotted:

pts2=Map[Drop[#,-1]&,pts3];

g2=Graphics[{
Map[Point[#]&,pts2]
}];

Show[g1,g2]

Oh, this is cool, it can be made into an animation! I learned this from Wolfram Research’s Theo Gray: http://www.youtube.com/watch?v=hhCpcw5Q9Zw

Each point is a “Locator” in the Manipulate[ ] function. The LocatorAutoCreate option permits additions and deletions from 3 to ∞. Use command-click or alt-click to create or delete new points. The colors pull-down is that ColorData[ ] option. 

A huge amount is packed into this!  But the result is spectacular. :-)

Manipulate[ListDensityPlot[MapIndexed[Append[#,#2[[1]]]&, points],PlotRange->{{0,10},{0,10}},
InterpolationOrder->0,ColorFunction->{ColorData[colors][#]&},
Mesh->All,
FrameTicks->False],
{{points,RandomReal[{0,10},{10,2}]},{0,0},{10,10},Locator,LocatorAutoCreate->{3,∞}},{colors, ColorData["Gradients"]}
]

Have fun! ---Joseph O’Rourke