Ordering Problem With VoronoiMesh

Imagine I want to model 2D cells moving in a tissue. Consider, as an example, Lloyd's relaxation algorithm. Using VoronoiMesh I'm able to mimic some cell movement. Now, if I colour one of the cells in a different manner, I can track it as the mesh moves.

However, this doesn't work as expected, since VoronoiMesh yields different cell ordering for different sets of points and therefore I get the wrong colouring in the wrong cell. Which was obtained by the following code, where the function VorR is inspired by this answer

rel = Function[{pts, zmp},
Block[{cells},
cells = MeshPrimitives[
VoronoiMesh[pts, {{-zmp, zmp}, {-zmp, zmp}}], "Faces"];
RegionCentroid /@
cells[[SparseArray[Outer[#2@#1 &, pts, RegionMember /@ cells, 1],
Automatic, False]["NonzeroPositions"][[All, 2]]]]]];
VorR = Function[{pt, s}, Module[{pts2, vor, vcells, mesh},
pts2 =
Flatten[Table[
TranslationTransform[{  2 s i, 2 s j}][pt], {i, -1, 1}, {j, -1,
1}], 2];
vor = VoronoiMesh[pts2];
vcells = Catenate[NearestMeshCells[{vor, 2}, #] & /@ pt];
mesh = MeshRegion[MeshCoordinates[vor], MeshCells[vor, vcells]]
]];

n = 20; T = 50;
pts = RandomReal[{-1, 1}, {n, 2}];
val = ReplacePart[ConstantArray[1, n], RandomInteger[n] -> 0.2];
vl = {pts};
vorl = {VorR[vl[], 1]};

For[i = 1, i <= T, i++,
sca = (i 0.01 + 1);
vl = Append[vl, rel[Last[vl], 1.05]];
vorl = Append[vorl, VorR[Last[vl]*sca, sca]];
];

colr = ConstantArray[Table[val[[j]], {j, n}], T];
cels = Table[
Table[Graphics[{RGBColor[0.5, 0.65, 0.5, colr[[j, p]]],
MeshPrimitives[vorl[[j]], 2][[p]]}], {p, n}], {j, T}];
ListAnimate[Table[Show[cels[[j]]], {j, T}]]

As it can be seen, the selected cell jumps between different places and this is due to the reordering of the cells (polygons) in MeshPrimitives[vorl[[i]],2], which occurs every time VoronoiMesh is applied. How can I solve this and get the correct ordering each time?

My attempt: Initially, I thought I could track the seeds updating positions (by Lloyd's algorithm) and associate each with a cell polygon, but that doesn't seem to work because the cells are not ordered by their generating seeds. Since movement is "relatively" slow, I then thought about tracking the position of the cell centroids in order to test, at each step, which cell permutation occurred and therefore correct the order at every step. And this seems to work, as seen by the following code (apologises for the "messy" code, there are definitely neater ways of writing this using Select for instance)

n = 20; T = 50;
pts = RandomReal[{-1, 1}, {n, 2}];
val = {ReplacePart[ConstantArray[1, n], RandomInteger[n] -> 0.2]};
vl = {pts};
vorl = {VorR[vl[], 1]};
vlc = {Map[RegionCentroid, MeshPrimitives[vorl[], 2]]};
val1 = {};
For[p = 1, p <= n, p++,
For[pi = 1, pi <= n, pi++,
If[RegionDimension[
RegionIntersection[MeshPrimitives[Last[vorl], 2][[p]],
Point[vlc[][[pi]]]]] =!= -Infinity,
val1 = Append[val1, val[[1, pi]]]
]]];

For[i = 1, i <= T, i++,

sca = (i 0.01 + 1);
vl = Append[vl, rel[Last[vl], 1.05]];
vorl = Append[vorl, VorR[Last[vl]*sca, sca]];
vlc = Append[vlc,
Map[RegionCentroid, MeshPrimitives[Last[vorl], 2]]];

nsol = {};
dsol = {};
For[p = 1, p <= n, p++,
If[RegionDimension[
RegionIntersection[MeshPrimitives[Last[vorl], 2][[p]],
Point[vlc[[i]][[p]]]]] =!= -Infinity,
nsol = Append[nsol, val[[i, p]]],
For[pi = 1, pi <= n, pi++,
If[RegionDimension[
RegionIntersection[MeshPrimitives[Last[vorl], 2][[p]],
Point[vlc[[i]][[pi]]]]] =!= -Infinity,
nsol = Append[nsol, val[[i, pi]]]
]]]];
val = Append[val, nsol]
];

colr = ReplacePart[val, 1 -> val1];
cels = Table[
Table[Graphics[{RGBColor[0.5, 0.65, 0.5, colr[[j, p]]],
MeshPrimitives[vorl[[j]], 2][[p]]}], {p, n}], {j, T}];
ListAnimate[Table[Show[cels[[j]]], {j, T}]] My previous code is even more clear if instead you take these values as

val = {Range[n]/n};

to get This, however, is not ideal, because it is based on the previous cell centroid, and whether it intersects the current cell in the current mesh. In other words, all I'm doing is the following: I test whether the previous centroid intersects the corresponding new cell (in the list ordering). If yes, I don't change the order, and if not I track which centroid is contained in this cell and reorder the values of each cell colour (val) accordingly. However, if the cell movement is too quick (per time iteration), the previous centroids are not guaranteed to intersect the new cells and order is compromised. So this is one problem

I wonder, then, if there is both a neater and more clever way of doing this, I wish Mathematica had this in mind, because even when defining, for example, adjacency matrices over a system of ODEs that describe some dynamics between cells, if the mesh is not static, then ordering becomes an recurrent issue, and even adding or removing cells (vertices) in the mesh (adjacency graph) requires careful management in order to get the correct ordering. I also noticed that reordering seems to take place whenever the number of neighbours of cells change (or degree of the corresponding graph vertex).

Any ideas? Sorry for the long post.

Note: the code doesn't always work due to Voronoi seeds being away from the cropping region, but it doesn't interfere with the main point of the question. Simply run it again. This was written in Mathematica 12.1.

• Code contains typo and not run at all in vv.12 .0-12.1.1.0. But as I understand you need some different code for modelling cells moving. May be not based on VoronoiMesh. Jun 30 '20 at 0:28
• There is no typo. Check again and check my note at the end. I need it to be based on VoronoiMesh, otherwise there is no point in making this question. Jun 30 '20 at 1:05
• Ok! I got it. Now I can imagine that you want "to model 2D cells moving in a tissue". But there is no any model of cells dynamic. So nobody can understand what the problem should be solved. Your solution is based on VoronoiMesh and Lloyd's relaxation algorithm. Do you try to improve this algorithm or just this code? Jun 30 '20 at 12:37
• Would you be interested in a solution to the cell dynamics problem that did not involve VoronoiMesh ? Jun 30 '20 at 12:58
• Could you generate a graph before/after you move the points by updating the GraphCoordinates option? As long as no new points are added or removed, this would keep a consistent ordering of the points, then you could apply the VoronoiMesh at every step, to get any of the properties of the cells that you might be interested in. Jun 30 '20 at 14:41

Update: I found that the slow part can be greatly improved by using the undocumented function RegionMeshMeshMemberCellIndex, as recommended in this question. The code is very similar to the previous version, but it runs much faster. For instance, the update of 100 points over 50 cycles of the Lloyd's algorithm takes about 15 seconds (as opposed to a couple minutes for ~16 points for 35 cycles, from the old version, running in a "normal" laptop).

(*How many cells?*)
n = 100;

(*Save consecutive {X,Y} coordinates here*)
spatialDomain = {-1, 1};
XYpositions = {RandomReal[spatialDomain, {n, 2}]};

(*How many time steps,for Lloyd's algorithm?*)
timeSteps = 50;

(*Ordering Array,this will be the correct indexing for the Voronoi \
cells*)
orderingArray = {};

i = 1;
While[i <= timeSteps,

(*Current XY positions,point coordinates*)
myPts = XYpositions[[-1]];

(*Current Mesh cells*)
currMesh = VoronoiMesh[myPts, {spatialDomain, spatialDomain}];
currMeshPrimitives = MeshPrimitives[currMesh, 2];

(*Correspondence Indexes between the current point orders and their \
mesh cell*)
Idx2 = #[] & /@ RegionMeshMeshMemberCellIndex[currMesh][myPts];

(*Append this to the Ordering array*)
AppendTo[orderingArray, Idx2];

(*Update the current XY points according to the mesh centroids*)
updateMeshCentroids =
RegionCentroid[#] & /@ currMeshPrimitives[[Idx2]];

(*Append the new XY points according to the correct order*)
AppendTo[XYpositions, updateMeshCentroids];
i++]

(*Choose some cell to "track"*)
trackThisCell = 20;
thisCellOverTime =
Table[orderingArray[[a]][[trackThisCell]], {a, 1,
Length[orderingArray]}];

And we get:

Manipulate[
VoronoiMesh[XYpositions[[a]], {spatialDomain, spatialDomain},
MeshCellLabel -> {2 -> "Index"},
MeshCellStyle -> {{2, _} -> LightBlue, {2, thisCellOverTime[[a]]} ->
LightGreen}], {a, 1, Length[XYpositions] - 1, 1}] Old version:

Here's a wildly inefficient way to do this, which nevertheless might be optimized/useful to you.

The main idea here is to identify if a given point is inside some cell in the Voronoi diagram before the transformation (in this way we ensure that no matter how "fast" the points move, we can "catch" them). This information is useful to know the identity of the cell after the transformation. To summarize the code below, we keep track of the correct index of every point to then map it to the corresponding cell in the Voronoi diagram(s).

We initialize some basic parameters and the arrays that will carry the useful information:

(*How many cells?*)
n = 16;

(*Save consecutive {X,Y} coordinates here*)
XYpositions = {RandomReal[{-1, 1}, {n, 2}]};

(*How many time steps, for Lloyd's algorithm?*)
timeSteps = 35;

(*Ordering Array, this will be the correct indexing for the Voronoi cells*)
orderingArray = {};

Now we run the process above described iteratively:

i = 1;
While[i <= timeSteps,

(*Current XY positions, point coordinates*)
myPts = XYpositions[[-1]];

(*Current Mesh cells*)
currMeshPrimitives =
MeshPrimitives[VoronoiMesh[myPts, {{-1, 1}, {-1, 1}}], 2];

(*Correspondence Indexes between the current point orders and their \
mesh cell*)

Idx = Flatten[
Table[Position[
RegionMember[#, myPts[[a]]] & /@ currMeshPrimitives, True], {a,
1, Length[myPts]}]];

(*Append this to the Ordering array*)
AppendTo[orderingArray, Idx];

(*Update the current XY points according to the mesh centroids*)
updateMeshCentroids =
RegionCentroid[#] & /@ currMeshPrimitives[[Idx]];

(*Append the new XY points according to the correct order*)
AppendTo[XYpositions, updateMeshCentroids];

i++] // AbsoluteTiming

So, in XYpositions we have the changes in the positions of the points, and in orderingArray we have the correct indexing of cells from this to the Voronoi cells.

Let's visualize one particular cell, say the 6th cell (note this is based on the identity of the points, not the current Voronoi cell label, which is the one that changes):

(*Choose some cell to "track"*)
trackThisCell = 6;
thisCellOverTime =
Table[orderingArray[[a]][[trackThisCell]], {a, 1,
Length[orderingArray]}];

To see that we are tracking a cell correctly, we can color it differently than the rest and see how it "moves". For comparison, I label the Voronoi cells with their "native" index, where you can see the problem of "inconsistent" labels across time (they change seemingly arbitrarily):

Table[VoronoiMesh[XYpositions[[a]], {{-1, 1}, {-1, 1}},
MeshCellLabel -> {2 -> "Index"},
MeshCellStyle -> {{2, _} -> LightBlue, {2, thisCellOverTime[[a]]} ->
LightGreen}], {a, 1, Length[XYpositions], 1}] I'm sure this code can be optimized, it runs slow mainly because of the way Idx is calculated. Although for a few dozen of cells is not bad. You might also need to implement a way to see if Lloyd's algorithm converges.

Sam,

I had this same issue a few years back and here is what I came up with. Let me just give you my bits and let you do the work of figuring out whether they work for your situation, but I believe they will.

Basically I adapted my functions from Quantum_Oli's answer at Find the nearest locations for multiple points

MatchTwoSetsOfPoints is the function you want. It is a wrapper for the more generalized MatchBallsToHoles which is a very nice and fast and non-statistical (which I believe means that it is comprehensive and perfect) routine for 'matching balls to holes', which is an assignment problem, and a special case of the 'minimum-cost flow problem'. The key functions are FindMinimumCostFlow and SourceTargetCostMatrix.

It also works for any dimensions of points.

Requires Mathematica v.10.2 for the FindMinimumCostFlow functions as used here. (for some reason AdjacencyGraph[costmatrix] doesn't work in 9.0).

There is a bug in FindMinimumCostFlow such that it sometimes takes days to evaluate ([CASE:4156292]), so I add a random factor to all elements with NudgeNonuniquePoints. Adding a random factor to ALL elements seems like overkill, it would be better to just add the random bits to the redundant points, but I don't bother.

SourceTargetCostMatrix is from Quantum_Oli; PositionsOfDuplicates is from Szabolcs; and GatherByList is from Woll on SE.

NudgeNonuniquePoints is by myself!

MatchTwoSetsOfPoints[balls_,holes_]:=("HolesOrdering"/.MatchBallsToHoles[balls,holes])/;Length[balls]==Length[holes]

PositionsOfDuplicates[list_List]:=DeleteCases[GatherByList[Range[Length[list]],list],{_}]

GatherByList[list_List,representatives_]:=Module[{funk},
funk/:Map[funk,_]:=representatives;GatherBy[list,funk]]

NudgeNonuniquePoints[ptsIn_,factor_:0.01]:=Module[{pts=ptsIn},
If[Length[pts]>Length[Union[pts]],
Map[Do[(pts[[elem]]=pts[[First[#]]]*(1+RandomReal[{-factor,factor},Dimensions[First[#]]])),{elem,Rest[#]}]&,PositionsOfDuplicates[pts]]];
pts]

SourceTargetCostMatrix[pointsA_,pointsB_]:=Module[{lA=Length[pointsA],lB=Length[pointsB]},ArrayFlatten@{{0,ConstantArray[1,{1,lA}],ConstantArray[0,{1,lB}],0},{ConstantArray[0,{lA,1}],ConstantArray[0,{lA,lA}],Outer[EuclideanDistance,pointsA,pointsB,1],ConstantArray[0,{lA,1}]},{ConstantArray[0,{lB,1}],ConstantArray[0,{lB,lA}],ConstantArray[0,{lB,lB}],ConstantArray[1,{lB,1}]},{0,ConstantArray[0,{1,lA}],ConstantArray[0,{1,lB}],0}}]

(*'FindMinimumCostFlow' requires mma10 for this use-case.*)
MatchBallsToHoles[ballsIn_,holesIn_]:=Module[{balls=ballsIn,holes=holesIn,nudge=0.01,costMatrix,assignments},
If[Length[holes]>Length[Union[holes]]||Length[balls]>Length[Union[balls]],Print["MatchBallsToHoles: WARNING: There were ",Length[balls]-Length[Union[balls]]," balls and ",Length[holes]-Length[Union[holes]]," holes that were in identical positions with other balls or holes that had to be perturbed by up to ",nudge*100," percent to avoid a bug in FindMinimumCostFlow."];];

(*'NudgeNonuniquePoints' is the 'Work-around' for when there are non-unique points that cause FindMinimumCostFlow to never converge:*)
balls=NudgeNonuniquePoints[balls,nudge];
holes=NudgeNonuniquePoints[holes,nudge];

costMatrix=SourceTargetCostMatrix[balls,holes];
assignments=Cases[FindMinimumCostFlow[costMatrix,1,Length[costMatrix],"EdgeList"],x_\[DirectedEdge]y_/;x!=1&&y!=Length[costMatrix]];

{"CostMatrix"->costMatrix,
"HolesOrdering"->assignments/.i_\[DirectedEdge]j_:>(j-Length[balls]-1),
"MatchedPoints"->assignments/.i_\[DirectedEdge]j_:>{balls[[i-1]],holes[[j-Length[balls]-1]]},
"NudgedBalls"->balls,"NudgedHoles"->holes}]
• .@berniethejet: if this is based on a nearest distance, wouldn't it be prone to mismatching points that happened to get closer to different point that was not really the point of interest but a neighbor point "passing by"? I originally thought of using a similar approach, but I think this would be a significant problem, especially with large step-displacements and many points. Jul 2 '20 at 10:40
• @TumbiSapichu yes, definitely. If we had some time-series history of movements then presumably we could set up some sort of predictions based on the those. But without that I don't know what other criteria could Pareto dominate all others. I believe this approach here works well with large numbers of points, but expect it won't work well with large-step displacements. See the link to the original post in my answer. Jul 2 '20 at 13:38
• Having some sort of 'time series' function that predicts the next point(s) based on the history of previous points will always add information to the matching, and it would be very interesting to combine the techniques in my answer with those of some timeseries models. However one will have to have some method of combining the predictions and in the end there will invariably be situations where one or the other (or both) models will make errors. In the end, the choice of which models to use will depend on the process one is trying to forecast. Jul 5 '20 at 17:42
• Thank you. This seems a very elegant answer, though I need some time to fully understand what is going on, I'm not familiar with some of the nomenclature. But I will try it. In the meantime I have found a way of doing it with the seeds (posted below), could be interesting for you as well. Jul 8 '20 at 10:54

I'm really happy with the solutions provided and they seem to do the trick.

Nonetheless, I'm sharing my solution with you. I managed to solve the ordering problem by tracking the generating seeds instead and defining a function per that translates the permutations occurring in the mesh cells every time there is an update of the seeds positions and number (and consequent Voronoi tessellation). With this I can update both the seed and val list order accordingly so that they match with the previous configuration.

The idea is the same as the intersecting centroids, but this time I guarantee each seed corresponds to the correct cell and therefore there is no risk of either overlapping or non-intersecting cells in the fast-moving mesh case. Here's the code

rel = Function[{pts, zmp},
Block[{cells},
cells = MeshPrimitives[
VoronoiMesh[pts, {{-zmp, zmp}, {-zmp, zmp}}], "Faces"];
RegionCentroid /@
cells[[SparseArray[Outer[#2@#1 &, pts, RegionMember /@ cells, 1],
Automatic, False]["NonzeroPositions"][[All, 2]]]]]];
VorR = Function[{pt, s}, Module[{pts2, vor, vcells, mesh},
pts2 =
Flatten[Table[
TranslationTransform[{  2 s i, 2 s j}][pt], {i, -1, 1}, {j, -1,
1}], 2];
vor = VoronoiMesh[pts2];
vcells = Catenate[NearestMeshCells[{vor, 2}, #] & /@ pt];
mesh = MeshRegion[MeshCoordinates[vor], MeshCells[vor, vcells]]]];
n = 20; T = 30; sca = 1; kl = {n};
pts = RandomReal[{-1, 1}, {n, 2}];
val0 = {Range[n]/n};
vl = {pts};
vorl = {VorR[vl[], 1]};
vll = {Table[
Select[Last[vl],
RegionDimension[
RegionIntersection[MeshPrimitives[Last[vorl], 2][[j]],
Point[#]]] =!= -Infinity &][], {j, Last[kl]}]};
per = Function[l,
Table[l[[j]], {j,
Table[Position[vl[], vll[[1, j]]][[1, 1]], {j, n}]}]];
val = {per[val0[]]};
For[i = 1, i <= T, i++,
vl = Append[vl, per[ rel[Last[vl], 1.05]]];
vorl = Append[vorl, VorR[sca Last[vl], sca]];
kl = Append[kl, Last[kl]];
vll = Append[vll,
Table[Select[Last[vl],
RegionDimension[
RegionIntersection[MeshPrimitives[Last[vorl], 2][[j]],
Point[#]]] =!= -Infinity &][], {j, Last[kl]}]];
per = Function[l,
Table[l[[j]], {j,
Table[Position[Last[vl], Last[vll][[j]]][[1, 1]], {j, n}]}]];
val = Append[val, per[Last[val]]]];
colr = val;
cels = Table[
Table[Graphics[{RGBColor[0.5, 0.65, 0.5, colr[[j, p]]],
MeshPrimitives[vorl[[j]], 2][[p]]}], {p, kl[[j]]}], {j, T}];
ListAnimate[Table[Show[cels[[j]]], {j, T}]] • I just updated my previous answer, using an undocumented function the performance increases many-fold. Just in case. Jul 12 '20 at 2:22