Ordering Problem With VoronoiMesh

Update: I found that the slow part can be greatly improved by using the undocumented function Region`Mesh`MeshMemberCellIndex, 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 = #[[2]] & /@ Region`Mesh`MeshMemberCellIndex[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}]