Delaunay triangulation, Voronoi diagrams and convex hulls in arbitrary dimension.
Delaunay triangulation, Voronoi diagrams and convex hulls.
Based on the qhull C library.
Consider this list of vertices (actually these are the vertices of a
polyhedron):
vertices = [[ -5, -5, 16 ] -- 0, [ -5, 8, 3 ] -- 1, [ 4, -1, 3 ] -- 2, [ 4, -5, 7 ] -- 3, [ 4, -1, -10 ] -- 4, [ 4, -5, -10 ] -- 5, [ -5, 8, -10 ] -- 6, [ -5, -5, -10 ] -- 7]
The delaunay function splits the polyhedron into simplices, the tiles of the
tesselation:
> import Delaunay> d <- delaunay vertices False False> _tiles dfromList[ ( 0, Tile{ _simplex =Simplex{ _points =fromList[ ( 2 , [ 4.0 , -1.0 , 3.0 ] ), ( 4 , [ 4.0 , -1.0 , -10.0 ] ), ( 5 , [ 4.0 , -5.0 , -10.0 ] ), ( 7 , [ -5.0 , -5.0 , -10.0 ] )], _circumcenter =[ -0.5000000000000009 , -3.0 , -3.499999999999999 ], _circumradius = 8.154753215150047, _volume = 78.0}, _neighborsIds = fromList [ 1 , 3 ], _facetsIds = fromList [ 0 , 1 , 2 , 3 ], _family = Nothing, _toporiented = False}), ( 1, Tile{ _simplex =......
The field _tiles is a map of Tile objects. The keys of the map are
the tiles identifiers. A Tile object has five fields:
_simplex, a Simplex object;
_neighborsIds, a set of tiles identifiers, the neighbors of the tile;
facetsIds, a set of facets identifiers, the facets of the tile;
family, two tiles of the same family share the same circumcenter;
toporiented, Boolean, whether the tile is top-oriented.
A Simplex object has four fields:
_points, the vertices of the simplex, actually a map of the vertices
identifiers to their coordinates
_circumcenter, the coordinates of the circumcenter of the simplex;
_circumradius, the circumradius;
_volume, the volume of the simplex (the area in dimension 2, the
length in dimension 1).
Another field of the output of delaunay is _tilefacets:
> _tilefacets dfromList[ ( 0, TileFacet{ _subsimplex =Simplex{ _points =fromList[ ( 4 , [ 4.0 , -1.0 , -10.0 ] ), ( 5 , [ 4.0 , -5.0 , -10.0 ] ), ( 7 , [ -5.0 , -5.0 , -10.0 ] )], _circumcenter = [ -0.5000000000000009 , -3.0 , -10.0 ], _circumradius = 4.924428900898053, _volume = 36.0}, _facetOf = fromList [ 0 ], _normal = [ 0.0 , 0.0 , -1.0 ], _offset = -10.0}), ( 1, TileFacet{ _subsimplex =......
This is a map of TileFacet objects. A tile facet is a subsimplex. The keys of
the map are the identifiers of the facets.
A TileFacet object has four fields: _subsimplex, a Simplex object,_facetOf, the identifiers of the tiles this facet belongs to (a set of one
or two integers), _normal, the normal of the facet, and offset, the offset
of the facet.
Finally, the output of delaunay has a _sites field, the vertices with
additional information:
> _sites dfromList[ ( 0, Site{ _point = [ -5.0 , -5.0 , 16.0 ], _neighsitesIds = fromList [ 1 , 3 , 7 ], _neighfacetsIds = fromList [ 15 , 16 , 17 ], _neightilesIds = fromList [ 5 ]}), ( 1, Site......
This is a map of Site objects. The keys of the map are the identifiers of
the vertices. A Site object has four fields:
_point, the coordinates of the vertex;
_neighsitesIds, the identifiers of the connected vertices;
_neighfacetsIds, a set of integers, the identifiers of the facets the
vertex belongs to;
_neightilesIds, the set of the identifiers of the tiles the vertex belongs
to.
The library allows to get the Voronoi diagram of a list of sites (vertices)
from the Delaunay tesselation. Here is a 3D example.
centricCuboctahedron :: [[Double]]centricCuboctahedron = [[i,j,0] | i <- [-1,1], j <- [-1,1]] ++[[i,0,j] | i <- [-1,1], j <- [-1,1]] ++[[0,i,j] | i <- [-1,1], j <- [-1,1]] ++[[0,0,0]]import Delaunayimport Voronoi3Dd <- delaunay centricCuboctahedron False Falsev = voronoi3 d
In some circumstances, one has to run the Delaunay tesselation including the
degenerate tiles in order to get the correct Voronoi diagram, that is to saydelaunay vertices False True.
The output of voronoi3 is a list of Voronoi cells given as pairs, each pair
consisting of a site and a list of edges.
This is the cell of the center [0, 0, 0]:
> last v( [ 0.0 , 0.0 , 0.0 ], [ Edge3 ( ( -0.5 , -0.5 , 0.5 ) , ( 0.0 , 0.0 , 1.0 ) ), Edge3 ( ( -0.5 , -0.5 , 0.5 ) , ( 0.0 , -1.0 , 0.0 ) ), Edge3 ( ( -0.5 , -0.5 , 0.5 ) , ( -1.0 , 0.0 , 0.0 ) ), Edge3 ( ( -0.5 , 0.5 , 0.5 ) , ( 0.0 , 0.0 , 1.0 ) ), Edge3 ( ( -0.5 , 0.5 , 0.5 ) , ( 0.0 , 1.0 , 0.0 ) ), Edge3 ( ( -0.5 , 0.5 , 0.5 ) , ( -1.0 , 0.0 , 0.0 ) ), Edge3 ( ( 0.5 , -0.5 , 0.5 ) , ( 0.0 , 0.0 , 1.0 ) ), Edge3 ( ( 0.5 , -0.5 , 0.5 ) , ( 0.0 , -1.0 , 0.0 ) ), Edge3 ( ( 0.5 , -0.5 , 0.5 ) , ( 1.0 , 0.0 , 0.0 ) ), Edge3 ( ( 0.5 , 0.5 , 0.5 ) , ( 0.0 , 0.0 , 1.0 ) ), Edge3 ( ( 0.5 , 0.5 , 0.5 ) , ( 0.0 , 1.0 , 0.0 ) ), Edge3 ( ( 0.5 , 0.5 , 0.5 ) , ( 1.0 , 0.0 , 0.0 ) ), Edge3 ( ( -0.5 , -0.5 , -0.5 ) , ( 0.0 , 0.0 , -1.0 ) ), Edge3 ( ( -0.5 , -0.5 , -0.5 ) , ( 0.0 , -1.0 , 0.0 ) ), Edge3 ( ( -0.5 , -0.5 , -0.5 ) , ( -1.0 , 0.0 , 0.0 ) ), Edge3 ( ( -0.5 , 0.5 , -0.5 ) , ( 0.0 , 0.0 , -1.0 ) ), Edge3 ( ( -0.5 , 0.5 , -0.5 ) , ( 0.0 , 1.0 , 0.0 ) ), Edge3 ( ( -0.5 , 0.5 , -0.5 ) , ( -1.0 , 0.0 , 0.0 ) ), Edge3 ( ( 0.5 , -0.5 , -0.5 ) , ( 0.0 , 0.0 , -1.0 ) ), Edge3 ( ( 0.5 , -0.5 , -0.5 ) , ( 0.0 , -1.0 , 0.0 ) ), Edge3 ( ( 0.5 , -0.5 , -0.5 ) , ( 1.0 , 0.0 , 0.0 ) ), Edge3 ( ( 0.5 , 0.5 , -0.5 ) , ( 0.0 , 0.0 , -1.0 ) ), Edge3 ( ( 0.5 , 0.5 , -0.5 ) , ( 0.0 , 1.0 , 0.0 ) ), Edge3 ( ( 0.5 , 0.5 , -0.5 ) , ( 1.0 , 0.0 , 0.0 ) )])
This is a bounded cell: it has finite edges only. The other ones are not
bounded, they have infinite edges:
> head v( [ -1.0 , -1.0 , 0.0 ], [ Edge3 ( ( -0.5 , -0.5 , 0.5 ) , ( 0.0 , -1.0 , 0.0 ) ), Edge3 ( ( -0.5 , -0.5 , 0.5 ) , ( -1.0 , 0.0 , 0.0 ) ), IEdge3( ( -0.5 , -0.5 , 0.5 ), ( -0.5773502691896258 , -0.5773502691896258 , 0.5773502691896258 )), Edge3 ( ( -0.5 , -0.5 , -0.5 ) , ( 0.0 , -1.0 , 0.0 ) ), Edge3 ( ( -0.5 , -0.5 , -0.5 ) , ( -1.0 , 0.0 , 0.0 ) ), IEdge3( ( -0.5 , -0.5 , -0.5 ), ( -0.5773502691896258 , -0.5773502691896258 , -0.5773502691896258 )), IEdge3 ( ( -1.0 , 0.0 , 0.0 ) , ( 1.0 , 0.0 , 0.0 ) ), IEdge3 ( ( 0.0 , -1.0 , 0.0 ) , ( 0.0 , -1.0 , 0.0 ) )])
The convexHull function of the ConvexHull module generates the convex hull
of a list of points.
import ConvexHullimport ConvexHull.Examples -- for the function randomInCubepoints <- randomInCube 100 -- 100 random points in a cubehull <- convexHull points False False Nothing
The vertices of the convex hull are stored in the field _hvertices:
> _hvertices hullfromList[ ( 3, Vertex{ _point =[ 0.7872072051657094 , 0.450772463858757 , 1.9900427529711773e-2 ], _neighfacets = fromList [ 42 , 43 , 47 , 48 ], _neighvertices = fromList [ 1 , 11 , 64 , 88 ], _neighridges = fromList [ 70 , 71 , 72 , 77 ]}), ( 6, Vertex......
The edges in the field _hedges:
> _hedges hullfromList[ ( Pair 14 70, ( [ 0.9215432980174852 , 0.8554065771602318 , 0.9842902519648512 ], [ 0.9497713758656887 , 0.998006476041318 , 0.7243639875028591 ])), ( Pair 84 99......
The facets in the field _hfacets:
> _hfacets hullfromList[ ( 0, Facet{ _fvertices =fromList[ ( 4, [ 1.5757133629105136e-3, 0.6442797662244039, 0.7058559215899725]), ( 67, [ 2.7500520534961326e-2, 0.37516259577251554, 0.7331611715042575]), ( 77, [ 3.46399386146774e-2, 5.575911794526589e-2, 0.46787034305814157])], _fridges =fromList[ ( 0, Ridge{ _rvertices =fromList[ ( 4, [ 1.5757133629105136e-3, 0.6442797662244039, 0.7058559215899725]), ( 77, [ 3.46399386146774e-2, 5.575911794526589e-2, 0.46787034305814157])], _ridgeOf = fromList [ 0 , 4 ]}), ( 1, Ridge{ _rvertices =fromList[ ( 4, [ 1.5757133629105136e-3, 0.6442797662244039, 0.7058559215899725]), ( 67, [ 2.7500520534961326e-2, 0.37516259577251554, 0.7331611715042575])], _ridgeOf = fromList [ 0 , 2 ]}), ( 2, Ridge{ _rvertices =fromList[ ( 67, [ 2.7500520534961326e-2, 0.37516259577251554, 0.7331611715042575]), ( 77, [ 3.46399386146774e-2, 5.575911794526589e-2, 0.46787034305814157])], _ridgeOf = fromList [ 0 , 1 ]})], _centroid =[ 2.1238724170849748e-2 , 0.3584004933140618 , 0.6356291453841239 ], _normal =[ -0.9930268604214181, -8.766369712550202e-2, 7.882087723357102e-2], _offset = 2.40848904384814e-3, _area = 4.0339144929987907e-2, _neighbors = fromList [ 1 , 2 , 4 ], _family = None, _fedges =fromList[ ( Pair 4 67, ( [ 1.5757133629105136e-3, 0.6442797662244039, 0.7058559215899725], [ 2.7500520534961326e-2, 0.37516259577251554, 0.7331611715042575])), ( Pair 67 77, ( [ 2.7500520534961326e-2, 0.37516259577251554, 0.7331611715042575], [ 3.46399386146774e-2, 5.575911794526589e-2, 0.46787034305814157])), ( Pair 4 77, ( [ 1.5757133629105136e-3, 0.6442797662244039, 0.7058559215899725], [ 3.46399386146774e-2, 5.575911794526589e-2, 0.46787034305814157]))]}), ( 1, Facet......
import HalfSpacesimport Data.Ratio ((%))x = newVar 1y = newVar 2z = newVar 3constraints =[ x .>= 0 -- shortcut for x .>=. cst 0, x .<= 3, y .>= 0, y .<=. cst 2 ^-^ (2%3)*^x, z .>= 0, z .<=. cst 6 ^-^ 2*^x ^-^ 3*^y ]
> hsintersections constraints False[ [ -1.1102230246251565e-16 , -1.1102230246251565e-16 , 6.0 ], [ 0.0 , 2.0 , 0.0 ], [ 0.0 , 0.0 , 0.0 ], [ 3.0 , 0.0 , 0.0 ] ]
The convex hull of a curve on the sphere:
The Voronoi cell of a point inside the Utah teapot:
The Voronoi diagram of a projection of the truncated tesseract:
The Voronoi diagram of a cube surrounded by three perpendicular circles:
