Approximating Voronoi diagram without any distance checks

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. • The Max[] in your code is effectively the metric being used in your fake diagram, so I'm unclear about your "without doing any distance checks" requirement. – J. M. will be back soon Oct 3 '17 at 9:55
• @J.M., um, thank you for pointing this out. I'd like to also point out that Max[] is only needed when two expanding areas contact each other. It is not in any way a distance check. I could replace it with picking a color at random from the ones present. I can redo my code without Max[] if it's not clear. – Yuriy S Oct 3 '17 at 9:58
• Part of the reason your code works, if I'm reading it right, is that you're working entirely in the first quadrant, which allows Max[] to act as a (Chebyshev) metric. Try doing it if you're not at the first quadrant (say, centered at the origin) so that you have negative coordinates. – J. M. will be back soon Oct 3 '17 at 10:02
• @J.M., please see the edit. I don't really understand what you mean since the algorithm works just as well without Max[] – Yuriy S Oct 3 '17 at 10:15

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[]}, {dims[]}];
initpts = Transpose[{RandomInteger[{1, dims[]}, p], RandomInteger[{1, dims[]}, p]}];
MapIndexed[
With[{pt = #1, val = First@#2},
nbrs = Select[(pt + # &) /@ {{0, 0}, {-1, 0}, {1, 0}, {0, -1}, {0, 1}},
1 <= #[] <= dims[] && 1 <= #[] <= dims[] &];
(initgrid[[Sequence @@ #]] = val) & /@ nbrs] &, initpts];
ArrayPlot[initgrid] Then define the update function:

update[array_] := Module[{m = Length[array], n = Length[array[]]},
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 - 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];

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[
Table[
updaterule2[
TakeDrop[
Flatten[#[[i - 1 ;; i + 1, j - 1 ;; j + 1, 1]]], {5}], #[[i,
j, 2]]],
{i, 2, dims[] - 1}, {j, 2, dims[] - 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.

• Thank you! This is exactly what I wanted – Yuriy S Oct 3 '17 at 23:15
• Though the "neighbourhoods" seem to stop being circular on a larger scale (see the top right angle of the last animation, where we can see an octagonal shape), so I suppose it wouldn't work in the limit? – Yuriy S Oct 3 '17 at 23:20
• @YuriyS Very true. And there's a discrepancy between the approximation and the overlay on the top edge which seems to be a direct result. I think it can probably be fixed by altering the rule, but I'll get back to you on that. – aardvark2012 Oct 3 '17 at 23:27
• It's also important to remember that we're giving up continuity and distance, and still getting a pretty good result. There are still some quirks to it, but these will in general be pretty minor unless you start with a large, sparsely seeded grid. For more densely seeded arrays, I think this algorithm has more issues with closely placed initial points than with "octagonal circles". – aardvark2012 Oct 4 '17 at 0:29
• @YuriyS No, not distance. Memory perhaps, which might be bending the rules a little. The ownership layer just records a cell's first nonzero state (ie, which region it belongs to), and then only the owner cells in the neighbourhood affect its subsequent dynamics. The dynamics are completely specified by the rules in Wolfram's book. Without the memory, those rules would just tear the regions apart. – aardvark2012 Oct 8 '17 at 13:18

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[]];
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[];
{x, y} = Round[Rest[pt]];
statetable[[x, y]] = frozen;
bandheap[] = 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= {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.) 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] (*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}, 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. 