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

5 Answers 5

21
$\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

$\endgroup$
20
+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

$\endgroup$
3
  • 2
    $\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, 2016 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, 2016 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, 2016 at 6:27
15
$\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

$\endgroup$
1
13
$\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

$\endgroup$
2
  • $\begingroup$ Wouldn't creating a NearestFunction be better? I still prefer RegionMemberFunction since it's listable. $\endgroup$
    – RunnyKine
    Aug 29, 2014 at 17:47
  • $\begingroup$ @RunnyKine, agreed... $\endgroup$
    – kglr
    Aug 29, 2014 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, #] &]]
$\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.