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?

@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. — RunnyKine, Feb 01 '13 at 23:43
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. — Szabolcs, Jun 18 '18 at 13:22

score 21 · Answer 1 · edited Dec 17 '16 at 03:43

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

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) http://math.lbl.gov/voro++/ — saturasl, Jul 15 '14 at 06:25

score 16 · Answer 2 · answered Feb 02 '13 at 03:14

16

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

answered Feb 02 '13 at 03:14

Szabolcs

234,956
30
623
1,263

Yes, TetGenLink is the way to go. – Feb 02 '13 at 03:17
1

@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. – Szabolcs Feb 02 '13 at 03:18
well TetGen will not give you the Voronoi directly. You will need to construct the Voronoi from the tetrahedralization. – Feb 02 '13 at 05:55
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. – Rainer Sep 22 '13 at 22:02
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 :) – ybeltukov Sep 21 '14 at 12:13

Greg Hurst · Answer 3 · 2018-11-17T17:19:52.743

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

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[]. — J. M.'s missing motivation, Feb 24 '19 at 14:51

J. M.'s missing motivation · Answer 4 · 2016-12-18T06:42:07.603

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

score 7 · Answer 5 · edited Apr 13 '17 at 12:56

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?

score 7 · Answer 6 · edited Mar 25 '17 at 17:00

7

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.

edited Mar 25 '17 at 17:00

J. M.'s missing motivation

124,525
11
401
574

answered Feb 02 '13 at 03:37

Tuku

486
2
4

Looks relevant: http://xlr8r.info/mPower/pages/boundedCellVoronoi.html – DavidC Feb 02 '13 at 03:43
3

This is quite a bare-bones answer. Please consider extending it, and adding an example if possible. – Mr.Wizard Feb 02 '13 at 05:53
@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. – Szabolcs Feb 02 '13 at 18:11
1

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? – Szabolcs Feb 02 '13 at 20:25
@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. – Mr.Wizard Feb 03 '13 at 04:55
@Mr.Wizard I was also thinking of answers like this where I practically gave the algorithm, but didn't spend the time to implement it. So it got zero upvotes. I just didn't want you (or others) to misunderstand this as complaining about the upvotes---I am truly not rep motivated any more. Also, I do think that this qhull/mPower answer is still too valuable to "hide" it as a comment. Most likely I could get mPower working again with an hour or two of work. – Szabolcs Feb 03 '13 at 23:01
@Mr.Wizard What would you do in the same position, i.e. if you knew about mPower, strongly suspected that it can be fixed up with a little bit of work, but you were not willing to spend the time on it at the moment? Would you post a comment or an answer with the link but no full solution? If I'm a googler coming upon this question, I'd much rather see it as an answer than burrowed comment. Maybe it should be a community wiki answer explaining the situation and encouraging anyone who finishes fixing mPower to edit it and post a full tutorial? But a new SE user wouldn't do that. – Szabolcs Feb 03 '13 at 23:04
@Szabolcs I'll have to think about your argument. For now notice that I did not delete this as "not an answer" so my position is not about answer versus comment, but rather quality of answer and voting practice. – Mr.Wizard Feb 04 '13 at 01:01
@Szabolcs qhull works for the Delaunay triangularization using the MathLink program qh-math, but the Voronoi diagram requires one to extract results from qhull differently and that's not implemented in qh-math. (It can calculate the diagram directly; just reading out the data isn't there yet.) Support can probably be added without great difficulty but it doesn't look entirely trivial either. – Oleksandr R. Feb 04 '13 at 02:48
@Tuku, thanks for uploading this, but I can't get it to work on Windows. It gives me those two warnings like you said. But uploading your sample code I get the error: Error: convexHull: qconvex failed: no output produced. – RunnyKine Feb 05 '13 at 18:12
Ok, so I got it to work on Windows. For those using Windows like myself, you don't have to change anything as described above. Just copy the files in the "bin" folder from qhull and place them in the "src" folder under the "qhull" folder you created in the UserBaseDirectory. – RunnyKine Feb 05 '13 at 18:28

score 6 · Answer 7 · edited Dec 18 '16 at 04:06

6

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

edited Dec 18 '16 at 04:06

J. M.'s missing motivation

124,525
11
401
574

answered Feb 02 '13 at 01:12

Joseph O'Rourke

4,731
25
42

1

Actually it does exist in Mathematica, see my answer. – Szabolcs Feb 02 '13 at 03:23
Say, we have a tetrahedralization of the points and the circumcenter to the tets, how does on construct the Voronoi cell of these tetrahedron? – Feb 02 '13 at 05:57
The Voronoi cells are bounded by the planar polygons formed by the centers of the circumspheres which include a given link in the tetrahedralization. – Xerxes Feb 02 '13 at 06:35
1

@Szabolcs: I am pleased to learn of TetGen! – Joseph O'Rourke Feb 02 '13 at 14:46
@Xerxes, and what do you do at the boundary, where there is no neighboring circumsphere? – Feb 02 '13 at 20:32
@JosephO'Rourke, could you explain a bit if there is anything that could be improved such that you would have found this earlier? – Feb 02 '13 at 20:57
1

@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. – Joseph O'Rourke Feb 03 '13 at 00:43
@JosephO'Rourke, I was trying to ask two things: First, what would have helped you such that you would have found out about the existence of TetGenLink earlier that you had and, secondly, if I am not mistaken a Delaunay tetrahedralization is the dual of the Voronoi partitioning. Can you point me to some literature where the transformation from the tetrahedralization to the dual is discussed. I hope this is clearer now. If not let me know. – Feb 04 '13 at 08:35
@reubenko: The literature on 2D Voronoi diagrams suffices, and is everywhere. Each Delaunay edge $e$ connects two sites $a$ and $b$. The plane that is the perpendicular bisector of $e$ includes a face of the Voronoi cells of $a$ and $b$. Intersecting the halfspaces for each face of one site yields the convex polyhedral cell for that site. – Joseph O'Rourke Feb 04 '13 at 12:28
@JosephO'Rourke, thanks; I'll see if I can write an implementation for this. If I succeed I'll post it here. – Feb 04 '13 at 20:43
@Joseph The link in this post is broken. If possible would you update it, please? – Mr.Wizard Sep 21 '14 at 05:22

score 3 · Answer 8 · answered Sep 18 '23 at 17:38

3

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

answered Sep 18 '23 at 17:38

eldo

67,911
5
60
168

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

8 Answers8

Linked

Related