How can I define a 3D version of the built-in VoronoiDiagram (VoronoiMesh in V10) function?

Question

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?

Answer 1

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]]

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]]

---

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]]

Answer 2

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).

Answer 3

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]

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

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

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]}]

Answer 4

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]]

Answer 5

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 &)
]

Is that roughly what you want?

Answer 6

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.

Answer 7

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:

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

Answer 8

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}}]