36
$\begingroup$

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

Mathematica graphics

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}]}}
|improve this question
$\endgroup$
  • 2
    $\begingroup$ Good question. I was sure that one of vor["Properties"] would hold the answer but no luck... $\endgroup$ – Simon Woods Aug 29 '14 at 18:46
  • $\begingroup$ @SimonWoods Yeah, I thought so too. I was disappointed to find out it didn't. $\endgroup$ – RunnyKine Aug 29 '14 at 18:49
  • $\begingroup$ @SimonWoods Maybe "PointInFaces" can give some information.But the order is strange identically. $\endgroup$ – yode May 26 '16 at 9:50
  • $\begingroup$ In this post,I want to use RelationGraph[RegionMember[#2, #1] &, pts, poly] to make it. $\endgroup$ – yode May 26 '16 at 9:52
17
$\begingroup$

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

enter image description here

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.

enter image description here

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], 
 PointSize@0.02, Red, MeshCells[vor, 0][[in]], Opacity@0.2, 
 Yellow, MeshCells[vor, 2]}], PointSize@0.02, Point@pts}

enter image description here

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]]] &, 
     conn[[pin]]@"AdjacencyLists"]] & @@@ cells;
Show@{g, Graphics@{Orange, PointSize@0.02, Point@DeleteCases[p1, {}]}}

enter image description here

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

enter image description here

|improve this answer
$\endgroup$
17
+50
$\begingroup$

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

Region`Mesh`MeshMemberCellIndex[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,

Region`Mesh`MeshMemberCellIndex[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, Region`Mesh`MeshMemberCellIndex[vor][pt]]}]

enter image description here

|improve this answer
$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – kirma Mar 10 '16 at 21:32
  • $\begingroup$ @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? $\endgroup$ – ilian Mar 10 '16 at 21:48
  • 3
    $\begingroup$ 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. $\endgroup$ – kirma Mar 11 '16 at 6:27
13
$\begingroup$

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]

Mathematica graphics

|improve this answer
$\endgroup$
12
$\begingroup$ 
c = PropertyValue[{vor, 2}, MeshCellCentroid];
cents = Join @@ (Nearest @ pts /@c);
cs2 = Thread[{Point/@cents,MeshPrimitives[vor, 2]}];
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]]

enter image description here

|improve this answer
$\endgroup$
  • $\begingroup$ Wouldn't creating a NearestFunction be better? I still prefer RegionMemberFunction since it's listable. $\endgroup$ – RunnyKine Aug 29 '14 at 17:47
  • $\begingroup$ @RunnyKine, agreed... $\endgroup$ – kglr Aug 29 '14 at 17:50
1
$\begingroup$

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, #] &]]
|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.