# Site - Cell Correspondence in Voronoi Diagram obtained via VoronoiMesh

Consider the following:

pts = RandomReal[4, {10, 2}];
vor = VoronoiMesh[pts, {{0, 4}, {0, 4}}];

Graphics[{GraphicsComplex[
MeshCoordinates[vor], {Thick, Blue, MeshCells[vor, 1],
PointSize[0.02], Red, MeshCells[vor, 0], Opacity[0.2], Yellow,
MeshCells[vor, 2]}], PointSize[0.02], Point[pts[[1]]], Green,
MeshPrimitives[vor, 2][[1]]}]


Looking at the above, we see that the Voronoi sites (points that generate the Voronoi diagram) do not correspond to the Voronoi cells (in the order) generated by VoronoiMesh. Since VoronoiMesh discards the points used in generating the diagram, there's no obvious way I can see to align the Voronoi sites with their corresponding cells. This feels like an oversight to me, since one will have to do some complex post-processing to obtain such correspondence.

My question is: is there an easy way to obtain a site - cell alignment of the Voronoi diagram obtained via VoronoiMesh other than using e.g. RegionMember which seems like an unnecessary overkill when this information should be available in the Mesh object. So e.g. with the points above, an output of the form shown below is desirable:

Transpose[{pts, MeshCells[vor, 2]}]

{{{3.15429516, 1.68154241}, Polygon[{20, 11, 6, 15}]},
{{0.925806386, 3.57266721}, Polygon[{12, 3, 1, 13}]},
{{2.12535296, 3.88158116}, Polygon[{14, 5, 8, 18}]},
{{3.52419343, 0.689832144}, Polygon[{18, 8, 7, 17}]},
{{0.590442758, 1.53821808}, Polygon[{10, 4, 9, 6, 11}]},
{{2.70632597, 3.81614775}, Polygon[{19, 13, 1, 5, 14}]},
{{2.23994147, 1.5141697}, Polygon[{15, 6, 9, 2, 16}]},
{{3.68657472, 2.93505855}, Polygon[{16, 2, 3, 12, 21}]},
{{0.12127097, 2.63380289}, Polygon[{17, 7, 4, 10, 22}]},
{{2.21266451, 2.69765978}, Polygon[{1, 3, 2, 9, 4, 7, 8, 5}]}}

• Good question. I was sure that one of vor["Properties"] would hold the answer but no luck... Commented Aug 29, 2014 at 18:46
• @SimonWoods Yeah, I thought so too. I was disappointed to find out it didn't. Commented Aug 29, 2014 at 18:49
• @SimonWoods Maybe "PointInFaces" can give some information.But the order is strange identically.
– yode
Commented May 26, 2016 at 9:50
• In this post,I want to use RelationGraph[RegionMember[#2, #1] &, pts, poly] to make it.
– yode
Commented May 26, 2016 at 9:52

There is also a currently undocumented internal function that may be useful.

RegionMeshMeshMemberCellIndex[mr] generates a function which can be applied to list of points, giving for each pt the index of the (first encountered) highest-dimensional cell of mr containing pt. For example,

RegionMeshMeshMemberCellIndex[vor][pts]

(* {{2, 2}, {2, 4}, {2, 5}, {2, 3}, {2, 9}, {2, 1}, {2, 10}, {2, 7}, {2, 8}, {2, 6}} *)

pt = pts[[1]];
Graphics[{GraphicsComplex[MeshCoordinates[vor],
{Thick, Blue, MeshCells[vor, 1],
PointSize[0.02], Red, MeshCells[vor, 0],
Opacity[0.2], Yellow, MeshCells[vor, 2]}],
PointSize[0.02], Magenta, Point[pt],
Opacity[0.2], Green, MeshPrimitives[vor, RegionMeshMeshMemberCellIndex[vor][pt]]}]


• Having a correct correspondence between points and Voronoi cells (as actually quite strongly implied in the documentation) would be even better. Any chance there would be a fix for this? :I Commented Mar 10, 2016 at 21:32
• @kirma I think it is a reasonable suggestion, but to my knowledge there is no guarantee of a particular ordering at present. Where/how exactly is it implied in the documentation? Commented Mar 10, 2016 at 21:48
• Hmmh. Maybe I interpreted the second bullet point on "Details and Options" section of VoronoiMesh a bit too strongly. Correct ordering would make life a lot easier on many applications of VoronoiMesh nonetheless, and occassional non-ordering is a surprising feature to say the least. Commented Mar 11, 2016 at 6:27

Of course it's not good that Mathematica forget initial points for Voronoi mesh. May be it is a bug. However one can easily recover all generating points directly from the mesh. It's interesting from theoretical point of view.

Let's consider one point of the Voronoi mesh

There are three pairs of equal angles $\alpha,\beta,\gamma$ around this point. Therefore $\alpha+\beta+\gamma=\pi$. Hence we know $$\gamma = \pi-(\alpha+\beta)$$ where $\alpha+\beta$ is an angle between known vectors ${\bf r}_1$ and ${\bf r}_2$. Then we can rotate ${\bf r}_3$ by angle $\gamma$ (to the left or the right) and obtain a ray to one of the generating points. Intersection of rays from different points will give the position of a generating point.

Here it is my realization of this method.

SeedRandom[0, Method -> {"MKL", Method -> {"Sobol", "Dimension" -> 2}}];
(* Voronoi cells looks better with sobol low-discrepancy random sequence *)
pts = RandomReal[4, {50, 2}];
boundary = {{0, 4}, {0, 4}};
vor = VoronoiMesh[pts, boundary];

q = MeshCoordinates[vor];
conn = # + #\[Transpose] &@
SparseArray[# -> ConstantArray[1, Length@#], {1, 1} Max@#] &[# & @@@MeshCells[vor, 1]];
in = Flatten@Position[q, {x_, y_} /;
boundary[[1, 1]] < x < boundary[[1, 2]] && boundary[[2, 1]] < y < boundary[[2, 2]]];
(* points inside boundaries *)

g = Graphics@{GraphicsComplex[
MeshCoordinates@vor, {Thick, Blue, MeshCells[vor, 1],
[email protected], Red, MeshCells[vor, 0][[in]], [email protected],
Yellow, MeshCells[vor, 2]}], [email protected], Point@pts}


cells = MeshCells[vor, 2];
δ = {{1, 0}, {0, 1}};
e = {{0, -1}, {1, 0}};
c = e\[TensorProduct]δ - δ\[TensorProduct]e;
p1 = If[Length@# < 2, {}, LeastSquares@##] & @@
Transpose@With[{pin = # ⋂ in, p = #},
With[{a = c.(#4 - #).(#3 - #).(#2 - #)}, {a, a.#}] & @@@
MapIndexed[q[[Join[pin[[#2]], # ⋂ p, #~Complement~p]]] &,
Show@{g, Graphics@{Orange, [email protected], Point@DeleteCases[p1, {}]}}


Inside the second With # is one of the not-boundary points in the Voronoi mesh. #2 and #3 is two neigbour points in the same polygon. #4 is the neigbour point in another polygon. The matrix c make proper rotation to obtain a vector a, which is perpendicular to a ray.

Unfortunately there some generating points near edge, which we can not find by ray intersection (there are less then two rays). I mark these remain points by {} and denote by rem. We will find them by reflection of calculated points (denoted by calc).

rem = Flatten@Position[p1, {}];
calc = Complement[Range@Length@cells, rem];
edges = Partition[#, 2, 1, 1]~Join~Partition[Reverse@#, 2, 1, 1] & @@@ cells[[calc]];
p1[[rem]] = Mean /@ DeleteCases[
If[p1[[calc]][[#]] != {},
ReflectionTransform[Cross[Subtract @@ q[[edges[[##]]]]],
q[[edges[[##, 1]]]]]@p1[[calc]][[#]], {}] & @@@
DeleteCases[Join @@ Position[edges, #] & /@ Partition[#, 2, 1, 1], {}] & @@@
MeshCells[vor, 2][[Flatten@Position[p1, {}]]], {}];
Max@Abs[Sort@pts - Sort@p1]


We recover all points with high precision! Of course, the order of points is the same as in MeshCells.

3.9968*10^-15

While I await other answers, here is the RegionMember approach I mentioned:

cellSite[p_, reg_] := With[{rm = RegionMember[reg]}, {Point@Flatten@Pick[p, rm[p]], reg}]


Then:

cs = cellSite[pts, #] & /@ MeshPrimitives[vor, 2];


Visualize:

GraphicsGrid[
Partition[Graphics[{GraphicsComplex[MeshCoordinates[vor], {Thick, Blue, MeshCells[vor, 1],
PointSize[0.02], Red, MeshCells[vor, 0], Opacity[0.2], Yellow, MeshCells[vor, 2]}],
PointSize[0.02], #1, Opacity[0.2], Green, #2, Red}] & @@@ cs, 5], ImageSize -> 700]


c = PropertyValue[{vor, 2}, MeshCellCentroid];
cents = Join @@ (Nearest @ pts /@c);
cs2 == cs
(* True *)

Grid[Partition[Graphics[{GraphicsComplex[MeshCoordinates[vor],
{Thick, Blue, MeshCells[vor, 1],
PointSize[0.02], Red, MeshCells[vor, 0], Opacity[0.2],
Yellow, MeshCells[vor, 2]}],
PointSize[0.02], #, Opacity[0.2],Green, #2, Red}] &@@@ cs2,5]]


• Wouldn't creating a NearestFunction be better? I still prefer RegionMemberFunction since it's listable. Commented Aug 29, 2014 at 17:47
• @RunnyKine, agreed...
– kglr
Commented Aug 29, 2014 at 17:50

Another simple solution in my this post but very slow

pts = RandomReal[4, {10, 2}];
vor = VoronoiMesh[pts, {{0, 4}, {0, 4}}];
poly = MeshPrimitives[vor, 2];

AbsoluteTiming[
youWant =
List @@@ EdgeList@Quiet[RelationGraph[RegionMember, poly, pts]];
RegionMember @@@ youWant]


{0.414993,{True,True,True,True,True,True,True,True,True,True}}

But I don't know why I get some error informations when I don't use Quiet.

Another more fast method can do this

youWant = Quiet[Gather[Join[pts, poly], # || RegionMember[#2, #] &]]