Message Boards Message Boards

Random Snowflake Generator Based on Cellular Automaton

Posted 11 years ago

Some time ago one of my friends asked me whether it is possible to design a cellular automaton which can generate realistic snowflakes. I recall my crystallography and thermodynamics knowledge and came up a very simple yet impressive model.

The Regular Triangular Lattice

First of all, we are trying to simulate snowflake, which is a kind of hexagonal crystal. So it should be best to construct our CA on a regular hexagonal grid, i.e. regular triangular lattice.

We all know CellularAutomaton inherently works on rectangle lattices ("4-lattice" for short), so how can we deduce a triangular lattice ("3-lattice" for short) on it? Well, the differences between rect-lattice and triangular one is just a geometric transformation.

To demonstrate that, have a look at the following 4-lattice, with a blue square highlighting the range-1 Moore neighborhood:

Clearly there is always a hexagon (the green area) in this kind of neighborhood.

So forming a regular 3-lattice is as straightforward as doing a simple affine transformation (basically a shearing and a scaling):

So to take advantage of all the power of CellularAutomaton, all we have to do, is to use a following special 6-neighborhood stencil on rectangle lattices, meanwhile our model can be discussed and constructed on regular triangular lattice convieniently:

And after the calculation, we can perform the affine transformation with following functions to get a nice hexagonal grid picture.

 vertexFunc =
         Compile[{{para, _Real, 1}},
             Module[{center, ratio}, center = para[[1 ;; 2]];
                 ratio = para[[3]];
                 {Re[#], Im[#]} + {{1, -(1/2)}, {0,
                                         Sqrt[3]/2}}.Reverse[{-1, 1} center + {3, 0}] & /@ (ratio 1/
                                 Sqrt[3] E^(I ?/6) E^(I Range[6] ?/3))],
            RuntimeAttributes -> {Listable}, Parallelization -> True,
            RuntimeOptions -> "Speed"

displayfunc[array_, ratio_] :=
            Polygon[vertexFunc[Append[#, ratio]] & /@ Position[array, 1]]},
        Background -> ColorData["DeepSeaColors"][0]]

The Model

To construct the crystallization model, let's consider one of the 6-neighborhood stencil, where each cell represents a minimal crystal unit:

A simple model will need only 2 states: 0 for "It's empty", 1 for "There is a crystal unit". So by considering all (except the 000000 one, because we are generating ONE snowflake thus don't want a crystall randomly arises from void) 6-bit non-negative numbers, we can have a finite set of possible arrangements of the neighborhood:

stateSet = Tuples[{0, 1}, 6] // Rest

However, from the viewpoint of physics, any two arrangements which can be transformed into each other with only rotation and reflection should be considered as the same arrangement in the sense of their physical effects on the central cell (i.e. cell2,2) are the same:

So we should gather stateSet with above equivalence class:

 gatherTestFunc = Function[lst, Union[Join[
                     RotateLeft[lst, # - 1] & /@ Flatten[Position[lst, 1]],
                     RotateLeft[Reverse[lst], # - 1] & /@
                         Flatten[Position[Reverse[lst], 1]]
 stateClsSet = Sort /@ Gather[stateSet, gatherTestFunc[#1] == gatherTestFunc[#2] &];
stateClsSetHomogeneous = ArrayPad[#, {{0, 12 - Length@#}, {0, 0}}] & /@ stateClsSet;

Which turned out to be 12 classes in total:

Now from the viewpoint of cellular automaton, we need to establish a set of rules on how should any 6-neighborhood arrangement, i.e. those 12 kinds of equivalence classes, determine the state of the central cell.

There are 4 kinds of possible transformations on cell2,2: 0 --> 1 is called frozen, 0 --> 0 is remaining empty, 1 --> 1 is remaining frozen, and 1 --> 0 is called melten. To make things more interesting and to explore more possibilities, we can introduce probability here, so certain arrangement will give certain probabilities corresponding to the 4 kinds of transformations. But notice that because of the unitarity of probability, we have Prob(frozen) + Prob(0->0) = 1 and Prob(melten) + Prob(1->1) = 1, so only 2 of the 4 probabilities are independent. In the following, we'll choose Prob(frozen) and Prob(melten), and denote them as pFrozen and pMelten.

Back to physics / thermodynamics, those 24 probabilities, pFrozen and pMelten, can of corse be determined by serious physical models, or they can be chosen randomly just for fun. For example, an intuitive (and naive) idea would be to believe an empty cell nearby a sharp pointed end or with abundant moisture source will have a high pFrozen. (People who are interested in the serious physical models should not miss the Gravner-Griffeath Snowfakes model.)

Now we have the grid, the stencil, the neighborhood arrangement set and the transfer probabilities, we're offically ready to construct our cellular automaton rules.
Following the above discussion, the construction is straightforward. There are only two points which need to pay attention to. One is to keep in mind that the rule function is applied on the 3x3 stencil, so even cell1,1 and cell3,3 has nothing to do with our model, don't forget handling them. The second is to use a SeedRandom function to make sure same arrangement gives same result in same time step, otherwise the 6-fold rotational symmetry and 3 axes of reflection symmetry will both break!

 ruleFunc = With[{
                 stateClsSetHomogeneous = stateClsSetHomogeneous,
                 seedStore = RandomInteger[{0, 1000}, 1000],
                 pFreeze = {1,   0,     0.6,   0,     0.3,   0.15,   0,     0.2,   0,     0.2,   0,     0.8},
                 pMelt   = {0,   0.7,   0.5,   0.7,   0.7,   0.5,    0.3,   0.5,   0.3,   0.2,   0.1,   0  }
             Compile[{{neighborarry, _Integer, 2}, {step, _Integer}},
                Module[{cv, neighborlst, cls, rand},
                    cv = neighborarry[[2, 2]];
                    neighborlst = {#[[1, 2]], #[[1, 3]], #[[2, 3]], #[[3, 2]], #[[3,
                                        1]], #[[2, 1]]} &[neighborarry];
                    If[Total[neighborlst] == 0, cv,
                        cls = Position[stateClsSetHomogeneous, neighborlst][[1, 1]];
                        SeedRandom[seedStore[[step + 1]]];
                        rand = RandomReal[];
                        Boole@If[cv == 0, rand < pFreeze[[cls]], rand > pMelt[[cls]]]
                (*CompilationTarget -> "C",*)
                RuntimeAttributes -> {Listable}, Parallelization -> True,
                RuntimeOptions -> "Speed"

(Note: re-compile the rule function ruleFunc will give a different set of seedStore thus a different growth path.)

Now everything is ready, let's grow a snowflake from the beginning! emoticon

 dataSet = Module[{
                     initM = {{
                                     {0, 0, 0},
                                     {0, 1, 0},
                                     {0, 0, 0}
                                 }, 0},
                     rspec = {1, 1},
                    tmin = 0, tmax = 100, dt = 1},
                rule = {ruleFunc, {}, rspec};
                CellularAutomaton[rule, initM, {{tmin, tmax, dt}}]
                ]; // AbsoluteTiming

    Rotate[displayfunc[dataSet[[k]], .8], 90 °],
    {k, 1, Length[dataSet], 1},
    AnimationDirection -> ForwardBackward,
    AnimationRunning -> False, DisplayAllSteps -> True

Possible Improvements

We used a SeedRandom function in our CA rule function to force the 6-fold rotational symmetry and 3 axes of reflection symmetry, and performed the CA calculation on all cells. However, this so called D6 symmetry can (and should) be integrated into our model, which will saving 11/12 of the calculation. Also, the randomness of the growth path comes from seedStore, so to generate a new growth path, we have to re-compile the rule function. But with a improved model as described above, this constraint will no longer exist.

Open question 

Can we construct a well-organized structure (like the crystals) from a cellular automaton defined on an irregular grid? While I believe the answer is yes, the next question would be how?
POSTED BY: Silvia Hao
7 Replies
Posted 11 years ago
check out chpt. 6 in "Modeing Nature with Cellular Automata using Mathematica". also see
POSTED BY: Richard Gaylord
Nice!  Here's a related post I wrote based on Matthew Szudzik's generator:

To Chris:

It's a very interesting post, thanks!

Last night I was thinking, currently the main difficulties for running large scale CA would be the speed. Maybe a function written in Assembly or running on GPU will be great.

PS. I remember your speaking on twisted architecture on a (maybe 2009) conference in China, it's very impressive emoticon
POSTED BY: Silvia Hao
Silvia, this post is wonderful, not only idea- but also presentation-wise. Thanks a lot for doing this! I am very curious if your model is close to those mentioned in NKS Book on page 371 which I think was implemented by Ed Pegg in: Snowflake-Like Patterns. But there we have a range of rules, - interestingly I also found a demo done by Matthew Szudzik that uses a specific numerically given rule 10926: Snowflake Growth. I bet all of these use mapping of rectangular to hexagonal grid because they use built in CA function.

POSTED BY: Sam Carrettie
To Simon:

Thank you! And you're absolutely right! emoticon

To Sam:

Thank you for your approving! I did this work about 2 years ago, back to when I haven't known the works on the NKS book or the Demonstration site. But I do have read them before I composed this post. They're fascinating. The rules they used are much concise than my model. I think I should take time read through the NKS book!

PS. The rectangular to hexagonal mapping is not essential, the built-in CA function can be used on any irregular grid (i.e. a general Graph). Actually I think I would like to write something about this topic.
POSTED BY: Silvia Hao
Posted 11 years ago
Cool idea!

You can speed up ruleFunc quite a bit by precomputing all random values:
 (* seedStore = RandomInteger[{0, 1000}, 1000], *)
 randStore = RandomReal[1, 1000],
  SeedRandom[seedStore[[step + 1]]];
  rand = RandomReal[];
 rand = randStore[[step+1]];                  
POSTED BY: Simon Schmidt
And by adjusting the pFrozen and pMelten, we're effectively adjusting the model, which then will generate different styles of shape (e.g. more branchy or more solid). Here are some examples obtained from different parameters:

POSTED BY: Silvia Hao
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract