Approximating Voronoi diagram without any distance checks

Question

When reading the Wikipedia article on Voronoi diagrams I was fascinated by this animation:

https://en.wikipedia.org/wiki/Voronoi_diagram#/media/File:Voronoi_growth_euclidean.gif

I noticed that to create a Voronoi diagram we just need to pick a locally interacting system with suitable properties such that a ball expands with constant speed around each of the initial points, but when the two balls meet they "annihilate" each other and stop expanding at that point.

However, I haven't seen any implementation of Voronoi diagrams without at least some distance checks. The dirtiest one just checks distance for every point to each of the initial points and picks the color of the nearest one.

As far as I can see even for approximate Voronoi diagrams on a grid distance checks are involved. See this wonderful question for example: Animating a Voronoi Diagram.

---

I have tried a quick implementation using cellular automaton on a square grid, however this way I don't get an Euclidean distance diagram, since each point expands as a square.

The question is: can we make a cellular automata or another kind of algorithm to create an approximation to Euclidean Voronoi diagram for a random set of points without doing any distance checks?

(Ideally, the approximation should approach the exact diagram for the infinite grid).

---

Here's the code I used:

Xm = 100;
Ym = 100;
Tm = 50;
Sm = 19;
F0 = Table[0, {j, 1, Xm}, {k, 1, Ym}];
Sx = Table[RandomInteger[{1, Xm}], {s, 1, Sm}];
Sy = Table[RandomInteger[{1, Ym}], {s, 1, Sm}];
Cs = Table[RandomInteger[{1, Sm^2}], {s, 1, Sm}];
Do[F0[[Sx[[s]], Sy[[s]]]] = -1, {s, 1, Sm}];
Do[l = If[Sx[[s]] == 1, Xm, Sx[[s]] - 1];
  r = If[Sx[[s]] == Xm, 1, Sx[[s]] + 1];
  t = If[Sy[[s]] == 1, Ym, Sy[[s]] - 1];
  b = If[Sy[[s]] == Ym, 1, Sy[[s]] + 1];
  F0[[l, Sy[[s]]]] = 
   F0[[r, Sy[[s]]]] = 
    F0[[Sx[[s]], t]] = F0[[Sx[[s]], b]] = Cs[[s]], {s, 1, Sm}];
F1 = F0;
ArrayPlot[F1, ColorFunction -> "ThermometerColors"]
Do[F0 = F1;
  Do[l = If[j == 1, Xm, j - 1];
   r = If[j == Xm, 1, j + 1];
   t = If[k == 1, Ym, k - 1];
   b = If[k == Ym, 1, k + 1];
   If[F0[[j, k]] == 0, 
    F1[[j, k]] = 
     Max[F0[[l, k]], F0[[r, k]], F0[[j, t]], F0[[j, b]]]], {j, 1, Xm}, {k, 1, Ym}], 
{n, 1, Tm}];
ArrayPlot[F1, ColorFunction -> "ThermometerColors"]

---

As a reaction to the comments: here is the relevant part of the code without using Max[]:

Do[l = If[j == 1, Xm, j - 1];
 r = If[j == Xm, 1, j + 1];
 t = If[k == 1, Ym, k - 1];
 b = If[k == Ym, 1, k + 1];
 If[F0[[j, k]] == 0 && 
   F0[[l, k]]^2 + F0[[r, k]]^2 + F0[[j, t]]^2 + F0[[j, b]]^2 > 0, 
  F1[[j, k]] = 
   RandomChoice[
    Cases[{F0[[l, k]], F0[[r, k]], F0[[j, t]], F0[[j, b]]}, 
     Except[0 | 0.]]]], {j, 1, Xm}, {k, 1, Ym}]

The only difference: we get a fuzzy boundary because of the uncertainty for the interacting colors. However, in the limit of grid size approaching infinity the fuzziness would disappear. No distance checks are involved.

---

I have redone the algorithm with a quasi hexagonal neighborhood (still using a square grid), but hexagons are not circles either, so I guess the main question is: how to make a cellular automaton which expands roughly as a circle (as the grid size and time tend to infinity)? If I can manage that, then I get an appoximation to a real Voronoi diagram.

---

Update

Using a method suggested by @aardvark2012 and some other tricks, I have been able to create expanding dodecagons (12-gons) which are very close to circles.

I need to iron my code a little before posting it here, but the basic idea is: I change the rule each 3 steps, so the direction of growth changes as well, on average I get this:

---

I haven't been able to draw the real Voronoi mesh over the diagram, so here's a crude juxtaposition (of the mesh created by Mathematica's procedure and the result of my code) made by hand. The disrepancy at the boundary is caused by the cyclic boundary conditions I use, not by the method itself.

Answer 1

Here's something that works (I think). The "trick" to getting quasi-circular neighbourhoods is to turn the initial points into "+"-like crosses, and then apply a "if three neighbouring cells have value x, set value to x" rule to the empty cells. Obviously you couldn't start from single cells with that rule, but once you get going it produces pretty good circles (until they start getting big... see comments):

Set up a grid with dimensions dims and p seed points. There is a slight problem with this initialisation in that if two seed points are located next to each other then their initial crosses will overlap, and that will mess up the process. There may be ways to mitigate this, but ultimately it's probably the price you pay for discretising.

dims = {100, 100}; p = 10;
initgrid = Table[0, {dims[[1]]}, {dims[[2]]}];
initpts = Transpose[{RandomInteger[{1, dims[[1]]}, p], RandomInteger[{1, dims[[2]]}, p]}];
MapIndexed[
  With[{pt = #1, val = First@#2},
    nbrs = Select[(pt + # &) /@ {{0, 0}, {-1, 0}, {1, 0}, {0, -1}, {0, 1}}, 
      1 <= #[[1]] <= dims[[1]] && 1 <= #[[2]] <= dims[[2]] &];
    (initgrid[[Sequence @@ #]] = val) & /@ nbrs] &, initpts];
ArrayPlot[initgrid]

Then define the update function:

update[array_] := Module[{m = Length[array], n = Length[array[[1]]]},
  ArrayPad[
   Table[
    If[array[[i, j]] == 0,
     With[{maxtally = 
        MaximalBy[
         Tally[Flatten@array[[i - 1 ;; i + 1, j - 1 ;; j + 1]] /. {0 -> Nothing}], 
        Last]},
      If[Length[maxtally] == 1 && maxtally[[1, 2]] >= 3, maxtally[[1, 1]], 0]],
     array[[i, j]]],
    {i, 2, m - 1}, {j, 2, n - 1}],
   1]
  ]

Then apply it until the process terminates:

gridlist = FixedPointList[update, ArrayPad[initgrid, 1]];

Note that initgrid gets padded so that neighbourhoods don't hang off the edges of the array.

Then, with a bit of messing, we can overlay the real Voronoi diagram over the array.

voronoigrid = gridlist[[-1, 2 ;; -2, 2 ;; -2]];
meshplot = Rotate[Show[
    VoronoiMesh[initpts, {1, #} & /@ dims, PlotTheme -> "Lines", 
     MeshCellHighlight -> {{1, All} -> {Thick, Black}}], 
    Graphics[{PointSize[Large], Point[initpts]}], 
    ImageSize -> {Automatic, 400}], -\[Pi]/2];

Overlay[{ArrayPlot[voronoigrid, Frame -> False, 
  ColorFunction -> "Rainbow", ImageSize -> 400], meshplot}]
ListAnimate[ArrayPlot[#, ColorFunction -> "Rainbow"] & /@ gridlist]

Update 1: Better images from a larger array (same algorithm).

Update 2: Using a "circular" CA.

For completeness, I wanted to look at using Wolfram's "approximate circle" CA... Not that it gives a bad result (comparatively), but it's a bit "fragile" in that the different populations can interact in surprising ways. Moreover, while it does indeed look -- on average -- approximately circle-like after a large number of generations, it does go through a large number of very un-circular configurations, which can make interactions at borders interesting. And sometimes the perimeter just stops moving. And it's v e r y slow.

Anyway, here's the code. First, the update rule. I had to build in an "ownership layer" to the array to keep track of which population activated a cell first, and hence, who gets to keep it -- otherwise they'd start eating each other.

updaterule2[{{state_}, nbd_}, owner_] := Module[{counts = Counts@nbd},
  If[owner > 0,
   If[KeyExistsQ[counts, owner],
    {Piecewise[{{owner, counts[owner] == 3}, {0, counts[owner] >= 5}},
       state], owner},
    {state, owner}],
   With[{count3 = Select[KeyDropFrom[counts, 0], # == 3 &]},
    If[Length[count3] == 1, {First@Keys[count3], 
      First@Keys[count3]}, {0, 0}]]
   ]
  ]

Initializing the array is complicated by the fact that two points can't be closer than 8 cells (it messes up the starting configurations which just breaks the whole thing. Like I said; "fragile"). The smallest symmetric starting configurations that lead to circle-like behaviour are 7 by 7 circles, given by DiskMatrix[3] - DiskMatrix[2, 7].

dims = {200, 200}; p = 16; initrad = 3;
initpts = {RandomInteger[{4, #}] & /@ (dims - 3)};
TimeConstrained[
 While[Length[initpts] < p,
  newpt = RandomInteger[{4, #}] & /@ (dims - 3);
  If[Min[EuclideanDistance[newpt, #] & /@ initpts] > 2 initrad + 2, 
   AppendTo[initpts, newpt]]
  ], 1]
p = Length[initpts]

initgrid = ConstantArray[{0, 0}, dims];
initmask = DiskMatrix[initrad] - DiskMatrix[initrad - 1, 2 initrad + 1];

Do[With[{i = initpts[[k, 1]], j = initpts[[k, 2]]},
  initgrid[[i - initrad ;; i + initrad, j - initrad ;; j + initrad, 
    1]] = k initmask;
  initgrid[[i - initrad ;; i + initrad, j - initrad ;; j + initrad, 
    2]] = k DiskMatrix[initrad]],
 {k, p}]

Then churn through:

gridlist = NestWhileList[
   newgrid = ArrayPad[
      Table[
       updaterule2[
        TakeDrop[
         Flatten[#[[i - 1 ;; i + 1, j - 1 ;; j + 1, 1]]], {5}], #[[i, 
         j, 2]]],
       {i, 2, dims[[1]] - 1}, {j, 2, dims[[2]] - 1}],
      {{1, 1}, {1, 1}, {0, 0}}] &,
   initgrid, UnsameQ[#1, #3] &, 3, 250];

(That UnsameQ[#1, #3] was supposed to catch period 2 repeating configurations. Of course, when I ran it, it was stuck in a period 3 pattern for 1500 iterations...)

Here's the raw layers -- showing the ON/OFF states of the cells as the regions grow:

and here's the "ownership layer", with the corresponding Voronoi mesh (using the same code from before, but with voronoigrid = gridlist[[-1, 2 ;; -2, 2 ;; -2, 2]]):

Not bad. Not great, though.

Answer 2

Here is a method that combines usage of "fast marching" and Nearest from responses here. It is not as pretty as the one shown here but is probably more efficient computationally. The idea is to work outward from the initial points, finding distance to nearest initially closest point, and index of that point, as we proceed. We store an intermediate state at some given "rate" (that is, every n pixels).

I should remark that some of the code in the first link uses Compile with inlining of compiled functions in a way that is currently broken. This has been reported and fixed for a future release.

Anyway, here is the code. Slightly clumsy but it seems to work.

frozen = -1.;
unseen = 0.;
far = 30000.;
outofbounds = 100000.;
bigstate = 10000;
band = 0.;

findNearestIndexC =
  Compile[{{ll, _Real, 2}, {nonzeros, _Integer, 2}, {rate, _Integer}},
    Catch[Module[{hindex = 0, dist, j1, j2, nbrs, pt, x, y, x1, y1, 
      x2, y2, next, prev, done, cond = False, len, wid, hsize, 
      statetable, statetable2, heaptable, bandheap, nbr, iter = 0},
     len = Length[ll];
     wid = Length[ll[[1]]];
     hsize = len*wid;
     nf = Nearest[N[nonzeros]];
     nf2 = Nearest[N[nonzeros] -> Range[Length[nonzeros]]];
     statetable = Map[If[TrueQ[# == unseen], far, #] &, ll, {2}];
     statetable2 = ll;
     heaptable = Table[0, {len}, {wid}];
     bandheap = Table[{0., 0., 0.}, {hsize}];
     Do[If[statetable[[ii, jj]] >= 0., Continue[]];
      nbrs = {{ii, jj + 1}, {ii, jj - 1}, {ii - 1, jj}, {ii + 1, 
         jj}};
      Do[{x, y} = nbrs[[kk]];
       If[! (1 <= x <= len && 1 <= y <= wid && 
           statetable[[x, y]] == far), Continue[]];
       hindex++;
       statetable[[x, y]] = band;
       nbr = nf[{x, y}];
       dist = Norm[First[nbr] - {x, y}];
       statetable2[[x, y]] = First[nf2[{x, y}]];
       bandheap[[hindex]] = {dist, N[x], N[y]};
       j1 = hindex;
       While[(j2 = Floor[j1/2]) >= 1 && 
         bandheap[[j2, 1]] > bandheap[[j1, 1]],
        bandheap[[{j1, j2}]] = bandheap[[{j2, j1}]];
        {x1, y1} = Round[Rest[bandheap[[j1]]]];
        heaptable[[x1, y1]] = j1;
        j1 = j2;];
       heaptable[[x, y]] = j1, {kk, Length[nbrs]}], {ii, len}, {jj, 
       wid}];
     While[hindex > 0,
      pt = bandheap[[1]];
      {x, y} = Round[Rest[pt]];
      statetable[[x, y]] = frozen;
      bandheap[[1]] = bandheap[[hindex]];
      done = False;
      prev = 1; next = 1;
      {j1, j2} = 2*prev + {0, 1};
      While[j1 < hindex && ! done, If[j2 < hindex,
        If[TrueQ[bandheap[[j1, 1]] <= bandheap[[j2, 1]]], next = j1, 
         next = j2], next = j1];
       cond = bandheap[[prev, 1]] > bandheap[[next, 1]];
       If[TrueQ[cond], 
        bandheap[[{prev, next}]] = bandheap[[{next, prev}]];
        {x1, y1} = Round[Rest[bandheap[[prev]]]];
        heaptable[[x1, y1]] = prev;
        prev = next;
        {j1, j2} = 2*prev + {0, 1};, done = True];];
      {x1, y1} = Round[Rest[bandheap[[prev]]]];
      heaptable[[x1, y1]] = prev;
      nbrs = {{x, y + 1}, {x, y - 1}, {x - 1, y}, {x + 1, y}};
      Do[{x2, y2} = nbrs[[kk]];
       If[! (1 <= x2 <= len && 1 <= y2 <= wid && 
           statetable[[x2, y2]] == band), Continue[]];
       nbr = nf[{x2, y2}];
       dist = Norm[First[nbr] - {x2, y2}];
       statetable2[[x2, y2]] = First[nf2[{x2, y2}]];
       j1 = heaptable[[x2, y2]];
       bandheap[[j1]] = {dist, N[x2], N[y2]};
       While[(j2 = Floor[j1/2]) >= 1 && 
         bandheap[[j2, 1]] > bandheap[[j1, 1]], 
        bandheap[[{j1, j2}]] = bandheap[[{j2, j1}]];
        {x1, y1} = Round[Rest[bandheap[[j1]]]];
        heaptable[[x1, y1]] = j1;
        j1 = j2;];
       heaptable[[x2, y2]] = j1, {kk, Length[nbrs]}];
      hindex--;
      iter++;
      If[Mod[iter, rate] == 0, Sow[statetable2]];
      Do[{x2, y2} = nbrs[[kk]];
       If[! (0 < x2 <= len && 0 < y2 <= wid && 
           statetable[[x2, y2]] == far), Continue[]];
       hindex++;
       statetable[[x2, y2]] = band;
       nbr = nf[{x2, y2}];
       dist = Norm[First[nbr] - {x2, y2}];
       statetable2[[x2, y2]] = First[nf2[{x2, y2}]];
       bandheap[[hindex]] = {dist, N[x2], N[y2]};
       j1 = hindex;
       While[(j2 = Floor[j1/2]) >= 1 && 
         bandheap[[j2, 1]] > bandheap[[j1, 1]], 
        bandheap[[{j1, j2}]] = bandheap[[{j2, j1}]];
        {x1, y1} = Round[Rest[bandheap[[j1]]]];
        heaptable[[x1, y1]] = j1;
        j1 = j2;];
       heaptable[[x2, y2]] = j1, {kk, Length[nbrs]}];];
     Sow[statetable2]]], {{nf[__], _Real, 2}},
   CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
   CompilationTarget -> "C"];

Here is an example, on a 100 x 200 grid with 20 points chosen at random. Preprocessing includes extracting the positions of those points. Takes around a fraction of a second on my Linux desktop machine.

dims = 100;
np = 20;
posns = Transpose[{RandomChoice[Range[dims], np], 
    RandomChoice[Range[dims], np]}];
pts = Normal[SparseArray[posns -> -Range[np], {dims, dims}]];

AbsoluteTiming[
 voronoi = Reap[findNearestIndexC[pts, posns, 400]][[2, 1]];]

(* Out[24]= {0.14349, Null} *)

Much slower to do the animation though.

ListAnimate[
 Map[MatrixPlot[# - pts, ColorFunction -> "BrightBands", 
    ImageSize -> 400] &, voronoi]]

(I uploaded a gif that seemingly was animated, but I guess it lost its interest. I do not know how to fix that, might try to figure it out when I get a bit of time free.)

Answer 3

I'm posting the code for 12-gons as an answer because the OP is bloated enough as it is. First, disclaimer: I still don't know how to properly use Mathematica, and the code is not by any means optimal. This is also the reason why I can't understand half of the code aardvark2012 used in the answer, so can't improve my own.

Here's the code with comments for each part (I got x and y confused, so x is for vertical direction):

Xm = 300; (*Field size, actually vertical*)
Ym = 300; (*Field size, actually horizontal*)
Tm = 100; (*Number of time steps*)
Sm = 19; (*Number of seed points*)
F0 = Table[0, {j, 1, Xm}, {k, 1, Ym}]; (*Empty field*)
Sx = Table[RandomInteger[{1, Xm}], {s, 1, Sm}];  (*Seed coordinates*)
Sy = Table[RandomInteger[{1, Ym}], {s, 1, Sm}];  (*Seed coordinates*)
Do[If[(OddQ[Sx[[s]]] && OddQ[Sy[[s]]]) || (EvenQ[Sx[[s]]] && 
     EvenQ[Sy[[s]]]), Sy[[s]] = Sy[[s]], 
  If[Sy[[s]] < Ym, Sy[[s]] = Sy[[s]] + 1, Sy[[s]] = Sy[[s]] - 1]], {s,
   1, Sm}];   (*Making sure all the seed points are lying on a quasi \
hexagonal grid, see the diagram*)
Cs = Table[2 + s, {s, 1, Sm}]; (*Setting colors for every seed*)
Do[F0[[Sx[[s]], Sy[[s]]]] = -1, {s, 1, 
  Sm}];  (*Putting seed points on the field*)
Do[l = If[Sx[[s]] == 1, Xm, 
   Sx[[s]] - 
    1];  (*Setting up the neighbourhood and cyclic boundary \
conditions*)
 ll = If[Sx[[s]] == 1, Xm - 1, If[Sx[[s]] == 2, Xm, Sx[[s]] - 2]];
 r = If[Sx[[s]] == Xm, 1, Sx[[s]] + 1];
 rr = If[Sx[[s]] == Xm, 2, If[Sx[[s]] == Xm - 1, 1, Sx[[s]] + 2]];
 t = If[Sy[[s]] == 1, Ym, Sy[[s]] - 1];
 b = If[Sy[[s]] == Ym, 1, Sy[[s]] + 1];
 F0[[ll, Sy[[s]]]] = 
  F0[[rr, Sy[[s]]]] = 
   F0[[l, t]] = F0[[l, b]] = F0[[r, t]] = F0[[r, b]] = Cs[[s]], {s, 1,
   Sm}];  (*Making seeds into small hexagons*)
F1 = F0; (*Array to be updated*)
F2 = F0; (*Filled array for plotting*)
Do[F0 = F1; 
 F2 = F0;
 Do[l = If[j == 1, Xm, j - 1];
  r = If[j == Xm, 1, j + 1];
  t = If[k == 1, Ym, k - 1];
  b = If[k == Ym, 1, k + 1];
  If[F0[[j, k]] == 
     0 && (F0[[l, k]] == F0[[r, k]] == F0[[j, t]] == F0[[j, b]]), 
   F2[[j, k]] = 
    Max[F0[[l, k]], F0[[r, k]], F0[[j, t]], F0[[j, b]]]], {j, 1, 
   Xm}, {k, 1, 
   Ym}]; (*Making sure all the cells we are not using are filled if \
surrounded by the same color*)
 Do[l = If[j == 1, Xm, j - 1];
  ll = If[j == 1, Xm - 1, If[j == 2, Xm, j - 2]];
  r = If[j == Xm, 1, j + 1];
  rr = If[j == Xm, 2, If[j == Xm - 1, 1, j + 2]];
  t = If[k == 1, Ym, k - 1];
  b = If[k == Ym, 1, k + 1];
  Cond = If[
    Mod[n, 4] == 
     1, (F0[[ll, k]]^2 F0[[l, t]]^2 > 0 || 
      F0[[ll, k]]^2 F0[[l, b]]^2 > 0 || 
      F0[[rr, k]]^2 F0[[r, t]]^2 > 0 || 
      F0[[rr, k]]^2 F0[[r, b]]^2 > 0 || 
      F0[[l, t]]^2 F0[[r, t]]^2 > 0 || F0[[l, b]]^2 F0[[r, b]]^2 > 0),
     F0[[ll, k]]^2 + F0[[l, t]]^2 + F0[[rr, k]]^2 + F0[[r, t]]^2 + 
      F0[[l, b]]^2 + F0[[r, b]]^2 > 0];
  (*The above condition ensures that for 1 step only cells with at \
least 2 filled neighbours grow, 
  while for 3 steps all the cells with at least 1 filled neighbour \
grow*)
  If[F0[[j, k]] == 0 && Cond, 
   F1[[j, k]] = 
    RandomChoice[
     Cases[{F0[[ll, k]], F0[[rr, k]], F0[[l, t]], F0[[l, b]], 
       F0[[r, t]], F0[[r, b]]}, Except[0 | 0.]]]], {j, 1, Xm}, {k, 1, 
   Ym}], {n, 1, 
  Tm}]; (*And lastly, around 2 or more colors the empty cell randomly \
picks between them. I would use Max, but then some of the expanding \
areas start to eat each other*)
ArrayPlot[F2, ColorFunction -> "ThermometerColors", 
 AspectRatio -> 
  1/Sqrt[3]] (*Since our hexagonal field is deformed, the resulting \
picture has to be rescaled. So we actually approximating Voronoi \
diagram for {Sy[[s]],Sx[[s]]/Sqrt[3]}, remember that here x is \
vertical*)

---

One improvement I see is to use "true" hexagonal neighbourhood instead of the one I'm using which looks like this:

---

This is the result on a smaller scale, looks very rough around the edges, but this can be fixed with more thorough conditions for the meeting of 2 or more colors.

---

Finally made a true hexagonal grid version, though at large scales it looks similar.

• I did not create a wraparound boundary conditions this time, since it's too complicated and unnecessary in this case.

• By using simple enough conditions I made sure that a cell can only change color if it's surrounded by a single color (with some empty cells), as well as that none of the seed points are each other's neighbours.