Approximating Voronoi diagram without any distance checks
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.
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.)
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.
