41
$\begingroup$

Can anybody point me in a direction that will guide me to extend the VoronoiDiagram function in Mathematica to handle 3D (three dimensional) situations (i.e. points in 3D)? Any help will be greatly appreciated.

Update

The current VoronoiMesh function in V10 still does not compute 3D Voronoi diagrams, which is strange since both the DelaunayMesh and ConvexHullMesh functions work with 3D data sets. I guess the waiting continues.... For the Wolfram Research employees, any ideas on why this was omitted, and what is the timeframe for when this functionality will be included in Mathematica?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ see `ComputationalGeometryPackage $\endgroup$ Commented Feb 1, 2013 at 23:24
  • 4
    $\begingroup$ @kguler, the VoronoiDiagram function only takes 2D data points and outputs the Voronoi diagram of points in a plane. What I'm trying to create is an equivalent VoronoiDiagram that accepts 3D vectors. $\endgroup$ Commented Feb 1, 2013 at 23:43
  • $\begingroup$ I think that the v10+ MeshRegion does not support arbitrary 3D cells, which would be needed for a representation of a Voronoi tessellation in 3D. So the mesh region functionality may not help much :-( One possibility might be using qhull. It has a command line interface, which can be accessed with RunProcess. MATLAB has qhull built in. If you have MATLAB, you can call it with MATLink. $\endgroup$ Commented Jun 18, 2018 at 13:22

8 Answers 8

Reset to default
21
$\begingroup$

Update: Using Raster3D and a variation of func that returns 4-tuples

 data3C = RandomReal[1, {10, 6}]; func3C = Nearest[{#, #2, #3} -> {##4, .03} & @@@ data3C]; tbl3C = Table[ First[func3C[{x, y, z}]] // Quiet, {x, 0, 1, .01}, {y, 0, 1, .01}, {z, 0, 1, .01}];

Examples:

 Row[Labeled[Graphics3D[Raster3D[tbl3C, ColorFunction -> #, Method -> {"InterpolateValues" -> True}], Background -> Black, ImageSize -> 400, SphericalRegion -> True], #, Top] & /@ {Hue, RGBColor, (GrayLevel[#[[1]], .03] &)}, Spacer[5]]

Voronoi images

colorRules = Thread[# -> (ColorData[1, "ColorList"][[;; Length@#]])] &[ Flatten[tbl3C, 2] // DeleteDuplicates] /. RGBColor -> List; Row[Labeled[ Graphics3D[Raster3D[tbl3C /. colorRules, ColorFunction -> #, Method -> {"InterpolateValues" -> True}], Background -> Black, ImageSize -> 400, SphericalRegion -> True], #, Top] & /@ {(RGBColor[#[[1]], #[[2]], #[[3]], .01] &), (RGBColor[#[[1]], #[[2]], #[[3]], .03] &), (RGBColor[#[[1]], #[[2]], #[[3]], .05] &)}, Spacer[5]]

Voronoi images at various opacities

Using version-9 built-in Image3D with @Mr.Wizard's func:

data = RandomReal[1, {20, 4}]; func = Nearest[{#, #2, #3} -> #4 & @@@ data]; dta = Table[func[{x, y, z}] // Quiet, {x, 0, 1, .005}, {y, 0, 1, .005}, {z, 0, 1, .005}]; Grid[Partition[ Image3D[dta, ImageSize -> 350, ColorFunction -> #, Background -> Black] & /@ {"SunsetColorsOpacity", "RainbowOpacity", "WhiteBlackOpacity", (Append[Blend[{LightBlue, Blue, Yellow}, #], #] &)}, 2]]

Wizard's example

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ The free package Voro++ does all the tricks, but it is in C++. This package is used in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) math.lbl.gov/voro++ $\endgroup$ Commented Jul 15, 2014 at 6:25
16
$\begingroup$

Mathematica ships with a TetGen interface called TetGenLink. To learn how to use TetGenLink is a bit more work than using the usual Mathematica functions, so I am not going to post a complete solution now.

But the way to go is using TetGenLink. It can compute a Delaunay tetrahedral mesh, which is the dual of the Voronoi partitioning. TetGen can also compute Voronoi partitions, but I am not sure if this function is exposed in TetGenLink, you'd have to check.

The function to use is TetGenTetrahedralize and I think you need to read the TetGen docs to understand the second argument (those flags should be the same as the command line options to TetGen).

Improve this answer
$\endgroup$
5
  • $\begingroup$ Yes, TetGenLink is the way to go. $\endgroup$ Commented Feb 2, 2013 at 3:17
  • 1
    $\begingroup$ @ruebenko If you post an answer with the precise command, I'll delete this one. I'd have to keep reading the docs for a while to figure out the right syntax. $\endgroup$ Commented Feb 2, 2013 at 3:18
  • $\begingroup$ well TetGen will not give you the Voronoi directly. You will need to construct the Voronoi from the tetrahedralization. $\endgroup$ Commented Feb 2, 2013 at 5:55
  • $\begingroup$ Have you found any undocumented Mathematica feature of the TetGenLink regarding Voronoi partitioning yet? I was not able to find any command which does the trick. $\endgroup$ Commented Sep 22, 2013 at 22:02
  • $\begingroup$ Or one can use TetGen through LibraryLink. However when I have to obtain 3D Voronoi mesh I just write points to a text file p.node and run the sytem command tetgen -v p.node :) $\endgroup$ Commented Sep 21, 2014 at 12:13
12
$\begingroup$

Note that there's currently no way to represent a collection of 3D Voronoi mesh cells in a MeshRegion or BoundaryMeshRegion.

Here's a routine that takes the dual of the DelaunayMesh and returns an Association where the keys are the points and the values are their respective Voronoi cells.

pad[δ_][{min_, max_}] := {min, max} + δ(max-min){-1, 1} VoronoiCells[pts_] /; MatrixQ[pts, NumericQ] && 2 <= Last[Dimensions[pts]] <= 3 := Block[{bds, dm, conn, adj, lc, pc, cpts, hpts, hns, hp, vcells}, bds = pad[.1] /@ MinMax /@ Transpose[pts]; dm = DelaunayMesh[pts]; conn = dm["ConnectivityMatrix"[0, 1]]; adj = conn . Transpose[conn]; lc = conn["MatrixColumns"]; pc = adj["MatrixColumns"]; cpts = MeshCoordinates[dm]; vcells = Table[ hpts = PropertyValue[{dm, {1, lc[[i]]}}, MeshCellCentroid]; hns = Transpose[Transpose[cpts[[DeleteCases[pc[[i]], i]]]] - cpts[[i]]]; hp = MapThread[HalfSpace, {hns, hpts}]; BoundaryDiscretizeGraphics[#, PlotRange -> bds]& /@ hp, {i, MeshCellCount[dm, 0]} ]; AssociationThread[cpts, RegionIntersection @@@ vcells] ]

Example:

SeedRandom[10000]; pts = RandomReal[1, {10, 3}]; vc = VoronoiCells[pts]

enter image description here

Show[MapIndexed[ BoundaryMeshRegion[#, MeshCellStyle -> {1 -> {Black, Thick}, 2 -> {ColorData[112][First[#2]]}}] &, Values[vc] ]]

enter image description here

Show[ MapIndexed[ BoundaryMeshRegion[#, MeshCellStyle -> {1 -> Black, 2 -> {Opacity[0.5], ColorData[112][First[#2]]}}] &, Values[vc] ], Graphics3D[{PointSize[Large], Point[pts]}], Method -> {"RelieveDPZFighting" -> True} ]

enter image description here

Note that this works in 2D too:

SeedRandom[10000]; pts = RandomReal[1, {10, 2}]; vc = VoronoiCells[pts]; Show[MapIndexed[ BoundaryMeshRegion[#, MeshCellStyle -> {1 -> {Black, Thick}, 2 -> {ColorData[112][First[#2]]}}] &, Values[vc] ], Epilog -> {PointSize[Large], Point[pts]}]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ This looks very nice. I guess one can modify the line computing bds so that one can have custom bounds just like in the built-in VoronoiMesh[]. $\endgroup$ Commented Feb 24, 2019 at 14:51
8
$\begingroup$

Here is another ContourPlot3D[]-based method for generating an approximate Voronoi diagram. The idea is originally due to Quílez.

BlockRandom[SeedRandom["voronoi"]; pts = RandomReal[{-1, 1}, {32, 3}]]; nf = Nearest[pts]; (* gradient-normalized function *) vfun[x_?NumericQ, y_?NumericQ, z_?NumericQ] := With[{np = nf[{x, y, z}, 2]}, (EuclideanDistance[{x, y, z}, np[[2]]] - EuclideanDistance[{x, y, z}, np[[1]]])/ EuclideanDistance[Normalize[{x, y, z} - np[[2]]], Normalize[{x, y, z} - np[[1]]]]] With[{ε = 1/100}, (* tiny number, increase if you want to see gaps between cells *) ContourPlot3D[vfun[x, y, z] == ε, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, MaxRecursion -> 1, Mesh -> False, PlotPoints -> 35]]

approximate Voronoi diagram

Improve this answer
$\endgroup$
7
$\begingroup$

It would help if you gave an example of a 3D Voronoi diagram.

Perhaps you want something like this, using Nearest as I did here.

Warning: this is very slow and uses a lot of memory!

data = RandomReal[1, {20, 4}]; func = Nearest[{#, #2, #3} -> #4 & @@@ data]; ContourPlot3D[ func[{x, y, z}] // Quiet, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, ColorFunction -> (Hue@#3 &) ]

Mathematica graphics

Is that roughly what you want?

Improve this answer
$\endgroup$
7
$\begingroup$

mPower is what your are looking for, which interfaces with Qhull

I have used this package with Mathematica 7 and 8 regularly on both Mac and Windows. Let's do a testing installation with Mathematica 9 on Mac 10.8; you can just delete the folder directly after testing.

mPower 1.0 for Mathematica 6.0

Qhull 12.01 Mac Bindary

For the Windows binary, you can grab it directly from qhull.org.

Unzip both packages, create a folder qhull under mpower folder, then copy bin folder from unzipped Qhull package, and put it under the qhull folder you just created.

In the mPower folder, open mPower.m in Mathematica:

$QHULL::usage="$QHULL should contain the name (as a string) of the folder that contains the QHULL binaries on your system. $QHULL = ToFileName[{$UserBaseDirectory, "Applications", "qhull", "src"}];

modify it to 

$QHULL = ToFileName[{NotebookDirectory[], "qhull", "bin"}];

or you can use the absolute path directly

Then, create a new notebook file test.nb at the mPower folder

SetDirectory[NotebookDirectory[]] << "mPower.m"

You will get two warnings:

Warning: regtet binary not found. Expected location: /mPower-1.0-11-May-2008-1631/qhull/bin/ Warning: pwrvtx binary not found. Expected location: /mPower-1.0-11-May-2008-1631/qhull/bin/

Ignore them. You are ready to use the qhull interface now.

Run the following code for testing:

(* random 3D points *) points = RandomReal[1, {40, 3}]; (* construct 3D convexhull *) ch = convexHull[points, convexHullFormat -> {facetNormals -> True, facets -> True}]; (* generate facets for 3D convexhull *) facetIndices = facets /. ch; loop[alist_] := Append[alist, alist[[1]]]; loopedFacetIndices = loop /@ facetIndices; index2xyz[ijklist_] := points[[#]] & /@ ijklist; loopedFacetCoordinates = index2xyz /@ loopedFacetIndices; (* visulization 3D convexhull *) convexObject = Graphics3D[{PointSize[Large], Point[points], FaceForm[], EdgeForm[Blue], Polygon /@ loopedFacetCoordinates}, ImageSize -> 500]

Check the documentation on boundedCellVoronoi. You can run the example code by copying them into the test.nb you just created.

Improve this answer
$\endgroup$
11
  • $\begingroup$ Looks relevant: xlr8r.info/mPower/pages/boundedCellVoronoi.html $\endgroup$ Commented Feb 2, 2013 at 3:43
  • 3
    $\begingroup$ This is quite a bare-bones answer. Please consider extending it, and adding an example if possible. $\endgroup$ Commented Feb 2, 2013 at 5:53
  • $\begingroup$ @Mr.Wizard It's a link to a ready-made package that contains the functionality. The package seems to interface with a mature library, i.e. I expect it to be fast and accurate. I think this is a very useful answer and it leads the way to a solution that is much better than anything with Nearest. Lately I have the feeling that on this site we expect too much. Answers which don't give working code are not upvoted even if they give the key to the solution. The OP should really be able to handle it from here. If I were the OP, I'd focus on this answer first and accept as soon as I got it working. $\endgroup$ Commented Feb 2, 2013 at 18:11
  • 1
    $\begingroup$ I took a look at the package, and it seems it's not immediately trivial to set it up with the current version of qhull. So I agree with MrW, if anyone gets it working, please describe how to do it. @Tuku did you get it working? $\endgroup$ Commented Feb 2, 2013 at 20:25
  • $\begingroup$ @Szabolcs I don't discount that we ask a lot on this site, but I disagree that it is too much. Your reversal in the following comment is proof enough of this to me. Answers must be more than a signpost to the information needed to do it yourself, or half the Q's could be answered with a naked link to mathprogramming-intro.org. It is often very hard to determine what is or is not trivial until you attempt it yourself, and those who would vote on such as answer are negligent in their duty if they vote without the knowledge from doing so. And that would be expecting too much. $\endgroup$ Commented Feb 3, 2013 at 4:55
6
$\begingroup$

The extension to 3D is quite non-trivial, and although Wizard's contour plot is impressive, it is an approximation. What you want does not exist in Mathematica. I might suggest Manifold Lab, which has its own issues but some incredible capabilities as well:

Voronoi model

Source: ggaliens [link broken, cannot find new one].

Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ Actually it does exist in Mathematica, see my answer. $\endgroup$ Commented Feb 2, 2013 at 3:23
  • $\begingroup$ Say, we have a tetrahedralization of the points and the circumcenter to the tets, how does on construct the Voronoi cell of these tetrahedron? $\endgroup$ Commented Feb 2, 2013 at 5:57
  • $\begingroup$ The Voronoi cells are bounded by the planar polygons formed by the centers of the circumspheres which include a given link in the tetrahedralization. $\endgroup$ Commented Feb 2, 2013 at 6:35
  • 1
    $\begingroup$ @Szabolcs: I am pleased to learn of TetGen! $\endgroup$ Commented Feb 2, 2013 at 14:46
  • 1
    $\begingroup$ @reubenko: I don't understand your question; sorry. Let me say this. Moving from the Delaunay triangulation graph to the Voronoi diagram is a nontrivial computation. The edges of the DT graph give you bisector planes. To compute the Voronoi cells, you would have to compute the convex body delimited by these planes. It is theoretically simple but in practice not easy. This is why I am so impressed with the ManifoldLab capabilities. $\endgroup$ Commented Feb 3, 2013 at 0:43
4
$\begingroup$

VoronoiMesh was updated in V 13.1 and now also works in 3D

VoronoiMesh[RandomReal[{-1, 1}, {15, 3}], MeshCellLabel -> None, MeshCellHighlight -> {{1, All} -> {Thickness[0.002], Black}, {0, All} -> {PointSize[0.015], Red}}]

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.