Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford

This notebook demonstrates the convex hull, Voronoi diagram, and Delaunay tessellation of a point set.

Figures 13.12 and 13.13 were created with this notebook.

See Chapter 13: Facets for details.

Updated May 2008

(* -------- Execute this cell first ------- *)
(* All subsequent cells below must be executed in the order they appear *)

(* Load the necessary packages -- This was necessary in version 5.2 *)
(* <<"ComputationalGeometry`" *)

(* Output results to folder that this file lives in *)
(SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];) 

(* Here are some drawing primitive settings *)

(* Set radius for drawing data points *)
radius=0.15;
triwidth= 0.005;
tesswidth = 0.008;

(* Create a data set of {x,y} *)

Print["Input data points:"]
points={{4.4,14},{6.7,15.25},{6.9`,12.8},{2.1`,11.1`},{9.5`,14.9`},{13.2`,11.9`},{10.3`,12.3`},{6.8`,9.5`},{3.3`,7.7`},{0.6`,5.1`},{5.3`,2.4`},{8.45`,4.7`},{11.5`,9.6`},{13.8`,7.3`},{12.9`,3.1`},{11,1.1`},{6,1},{8,1},{10,10},{11,11}}

(* Find number of points in list *)
pointnum=Length[points];

Print["Number of input data points = ",pointnum];


(* Load the point set into a table of disks for drawing *)
(* This disk radius can be set in the previous cell *)
graphpts=Table[{Disk[points[[i]],radius]},{i,1,pointnum}];

(* Plot the data points *)
pointplot=Graphics[graphpts,AspectRatio->Automatic]

(* Export point set to postscript file *)
Export["points.eps",pointplot];

Play with the convex hull of the point set.

Imagine pegs are sticking out of the points above. If we put a rubber band around the outside of the pegs, then the polygon formed by the band forms the convex hull.

<<ComputationalGeometry`

(* Define the convex hull of the point set *)

Print["The points on the convex hull include "]
convexhull=ConvexHull[points]

(* Find number of points on convex hull *)
numch=Length[convexhull];

Print["Number of points on the convex hull = ",numch]

(* Notice -- to access a list element, double brackets are needed *)
chpolygon = Table[{points[[convexhull[[k]],1]], points[[convexhull[[k]],2]]}, {k,1,numch}]

(* Draw the convex hull *)
Graphics[{{Yellow, EdgeForm[Thick],Polygon[chpolygon]}, graphpts}]

Delaunay Triangulation

(* There are several ways to generate the triangulation and tessellation *)
(* Additional ones are commented out below. *)


(* Load the Voronoi diagram vertices into vorvert and the vertex adjacency list into vorval *)
(* {vorvert,vorval}=VoronoiDiagram[points]; *)

(* Load the Delaunay triangulation into delval *)
(* VoronoiDiagram[points,delval]; *)

(* Create the Delaunay triangulation of points *)
(* If you want to see the point labels (numbers), set flag to ture *)
Delaunay=PlanarGraphPlot[points, BaseStyle->{"FontSize"->30}, LabelPoints->False];

picdel=Show[{Delaunay,pointplot}]
Export["Delaunay.eps",picdel]

Voronoi Diagram

(* BaseStyle->{"FontSize"->1} *)

voronoi=DiagramPlot[points, BaseStyle->{"FontSize"->10}, LabelPoints->False];
picvor=Show[{voronoi,pointplot}]
Export["voronoi.eps",picvor];

Let's try a larger data set -- created using the random number generator

Question: What if two points are identical? What will Mathematica do? Give this a try.

pts2=Table[RandomReal[{0,1},2],{i,1,100}];
numpts2=Length[pts2];

(* Print out the point set *)
(*
Print["Point set with ",numpts2," pts2"];
Do[Print["Point ",i," ",pts2[[i,1]]," ",pts2[[i,2]]],{i,1,numpts2}]
*)

(* Create points as disks *)
diskpts2=Table[{Disk[pts2[[i]],0.01]},{i,1,numpts2}];

(* Plot the disk points *)
plotpts2=Graphics[diskpts2,AspectRatio->Automatic];



(* Define the convex hull of the point set *)

convexhull2=ConvexHull[pts2];

(* Find number of points on convex hull *)
numch2=Length[convexhull2];

(* Notice -- to access a list element, double brackets are needed *)
chpolygon2 = Table[{pts2[[convexhull2[[k]],1]], pts2[[convexhull2[[k]],2]]}, {k,1,numch2}]

(* Draw the convex hull *)
Graphics[{{Yellow, EdgeForm[Thick],Polygon[chpolygon2]}, diskpts2}] 


(* Delaunay triangulation *)
Delaunay=PlanarGraphPlot[pts2, BaseStyle->{"FontSize"->30}, LabelPoints->False];

picdel2=Show[{Delaunay,plotpts2}]

(* Voronoi Diagram *)
(* PlotRange Automatic clips outliers *)
voronoi2=DiagramPlot[pts2, BaseStyle->{"FontSize"->10}, LabelPoints->False];
picvor2=Show[{voronoi2,plotpts2}, AspectRatio->1, PlotRange->Automatic]