**DiscreteMath****`****ComputationalGeometry****`**

Computational geometry is the study of efficient algorithms for solving geometric problems. The "nearest neighbor" problem involves identifying one point out of a set of points that is nearest to the query point according to some measure of distance. The "nearest neighborhood" problem involves identifying the locus of points lying nearer to the query point than to any other point in the set. This package provides functions for solving these and related problems in the case of planar points and the Euclidean distance metric.

Computational geometry functions.

The convex hull of a set is the boundary of the smallest set containing . The Voronoi diagram of is the collection of "nearest neighborhoods" for each of the points in . For points in the plane, these neighborhoods are polygons. The Delaunay triangulation of is a triangulation of the points in such that no triangle contains a point of in its circumcircle. This is equivalent to connecting the points in

according to whether their neighborhood polygons share a common side.

This loads the package.
In[1]:= **<<DiscreteMath`ComputationalGeometry`**

Here is a list of points in the plane.
In[2]:= **data2D = {{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}};

This gives the indices of the points lying on the convex hull in counterclockwise order.
In[3]:= **convexhull = ConvexHull[data2D]**

Out[3]=

Duplicate points are ignored.
In[4]:= **ConvexHull[{{0,0}, {1,0}, {0,0}, {2,0},**

{1,1}}]

Out[4]=

This gives the counterclockwise vertex adjacency list for each point in the Delaunay triangulation. For example, the entry {1,{4,3,2}} indicates that the first point in data2D is connected in counterclockwise order to the fourth, third, and second points.
In[5]:= **(delval =**

DelaunayTriangulation[data2D]) // Short[#,2]&

Out[5]//Short=

While DelaunayTriangulation need only specify the connections between points, VoronoiDiagram must specify both a set of diagram vertices and the connections between those vertices. Another difference between the two functions is that while a triangulation consists of segments, a diagram consists of both segments and rays. For example, in the case of a Voronoi diagram, points in the interior of the convex hull will have nearest neighborhoods that are closed polygons, but the nearest neighborhoods of points on the convex hull will be open polygons.

These considerations make the output of VoronoiDiagram more complex than that of DelaunayTriangulation. The diagram is given as a list of diagram vertices followed by a diagram vertex adjacency list. The finite vertices of the diagram are listed first in the vertex list. The vertices lying at infinity have head Ray and are listed last.

This assigns the list of Voronoi diagram vertices to vorvert and the Voronoi diagram vertex adjacency list to vorval.
In[6]:= **{vorvert, vorval} = VoronoiDiagram[data2D];**

The first vertex in vorvert is a finite diagram vertex having coordinates {-0.0158537,8.44146}. The last vertex in vorvert is a point at infinity. This point is represented by a Ray object having origin {10.5172,3.46115} and containing {13.95,0.2}.
In[7]:= **vorvert // Short[#,2]&**

Out[7]=

Each entry in vorval gives the index of a point in data2D followed by a counterclockwise list of the Voronoi diagram vertices that comprise the point's nearest neighborhood polygon.
In[8]:= **vorval // Short**

Out[8]//Short=

Here is the Voronoi polygon vertex adjacency list for the first point in data2D.
In[9]:= **vorval[[1,2]]**

Out[9]=

This selects the coordinates of the polygon vertices from vorvert. The first two vertices have head List while the last two have head Ray. Thus, the Voronoi polygon associated with the first point in data2D is open and is defined by a segment and two rays.
In[10]:= **vorvert[[%]]**

Out[10]=

Computing the Voronoi diagram using the Delaunay triangulation and the convex hull.

This computes the Voronoi diagram of data2D more efficiently by making use of the Delaunay triangulation vertex adjacency list.
In[11]:= **VoronoiDiagram[data2D, delval];**

Here the Voronoi diagram is computed using both the Delaunay triangulation and the convex hull.
In[12]:= **VoronoiDiagram[data2D, delval,**

convexhull];

Computational geometry plotting functions.

The default of PlanarGraphPlot is a plot of the Delaunay triangulation of the points.
In[13]:= **PlanarGraphPlot[data2D,**

DefaultFont -> {"Courier", 8.}]

This plots the convex hull of the points.
In[14]:= **PlanarGraphPlot[data2D, convexhull]**

Here is an alternative triangulation.
In[15]:= **trival = Insert[Insert[Delete[delval,**

{{12, 2, 4}, {16, 2, 2}}], 15, {11, 2, 2}],

11, {15, 2, 3}];

This plots the triangulation of data2D given by trival.
In[16]:= **PlanarGraphPlot[data2D, trival,**

DefaultFont -> {"Courier", 8.}]

The default of DiagramPlot is a plot of the Voronoi diagram of the points.
In[17]:= **DiagramPlot[data2D]**

Here is an alternative set of diagram vertices.
In[18]:= **diagvert = ReplacePart[vorvert,**

{-6., 0.}, {27, 2}];

Here is an alternative diagram vertex adjacency list.
In[19]:= **diagval = Join[Drop[vorval, -8],**

{{9, {1, 6, 8, 3}}, {10, {2, 6, 1,

24, 27}}},

Drop[vorval, 10]];

This plots the diagram of data2D given by diagvert and diagval.
In[20]:= **DiagramPlot[data2D, diagvert, diagval]**

Here is a set of three-dimensional points having the same coordinates as data2D

.
In[21]:= **data3D = Map[Append[#,**

Sqrt[64-(#[[1]]-8)^2-(#[[2]]-8)^2]]&,

data2D];

The default of TriangularSurfacePlot is a plot of the coordinates according to the connectivity established by the Delaunay triangulation of the

coordinates.
In[22]:= **TriangularSurfacePlot[data3D]**

This plots the coordinates according to the connectivity established by the triangulation trival

.
In[23]:= **TriangularSurfacePlot[data3D, trival]**

Options for computational geometry functions.

This gives the minimum set of points needed to define the convex hull.
In[24]:= **ConvexHull[{{0,0}, {1,0}, {0,0}, {2,0},**

{1,1}}, AllPoints -> False]

Out[24]=

This returns both the Delaunay triangulation and the convex hull.
In[25]:= **DelaunayTriangulation[data2D,**

Hull -> True] // Shallow

Out[25]//Shallow=

Here is a set of random numbers uniformly distributed on [0,1]

[0,1].
In[26]:= **random = Table[{Random[], Random[]}, {40}];**

This computes the Voronoi diagram of random.
In[27]:= **{randvert, randval} =**

VoronoiDiagram[random];

The diagram plot is dominated by outlier vertices.
In[28]:= **DiagramPlot[random, randvert, randval,**

LabelPoints -> False]

TrimPoints->2 means that the diagram is plotted so that both the outermost ray and the outermost vertex are eliminated.
In[29]:= **DiagramPlot[random, randvert, randval,**

LabelPoints -> False, TrimPoints -> 2]

The TrimPoints option can be used to "magnify" the diagram until the points in random fill the plot.
In[30]:= **DiagramPlot[random, randvert, randval,**

LabelPoints -> False, TrimPoints -> 6]

Testing for a Delaunay triangulation.

delval is a Delaunay triangulation, so this returns True.
In[31]:= **DelaunayTriangulationQ[data2D, delval]**

Out[31]=

This returns False, since trival is not a Delaunay triangulation.
In[32]:= **DelaunayTriangulationQ[data2D, trival]**

DelaunayTriangulationQ::inval: Triangle {11, 15, 12} is not a valid Delaunay triangle.

Out[32]=

Computing the nearest neighbor.

This computes the point in data2D nearest to {7.92,8.92}.
In[33]:= **neighbor = NearestNeighbor[{7.92, 8.92},**

data2D]

Out[33]=

If the Voronoi diagram is known, the diagram vertex list and vertex adjacency list can be substituted for the point list for a faster NearestNeighbor calculation.
In[34]:= **NearestNeighbor[{7.92, 8.92},**

vorvert, vorval]

Out[34]=

This gives the coordinates of the Voronoi polygon containing {7.92,8.92}.
In[35]:= **vorvert[[vorval[[neighbor, 2]]]]**

Out[35]=

For each of the points in data2D, the nearest point in data2D is the point itself.
In[36]:= **NearestNeighbor[data2D, vorvert, vorval]**

Out[36]=

The first half of these points is known to derive from one distribution; the second half is known to derive from another.
In[37]:= **known = Join[{{5.84, 1.2}, {5.94, 10.99},**

{5.1, 2.82}, {5.8, 1.67}, {5.63, 10.}},

{{0.31, 5.11}, {7.73, 5.38},

{10.42, 5.89}, {6.1, 5.1}, {6.92, 5.63}}];

Each of these points is believed to derive from one of the two distributions.
In[38]:= **unknown = {{5.56, 7.48}, {5.1, 1.67},**

{5.17, 4.89}, {0.3, 5.27},

{6.74, 5.73}, {5.09, 9.07}};

This classifies the points in unknown according to the classifications of their nearest neighbors in known.
In[39]:= **If[# > 5, 2, 1]& /@**

NearestNeighbor[unknown, known]

Out[39]=