Higher-dimensional generalizations [of substitution systems]

The state of a d-dimensional substitution system can be represented by a nested list of depth d. The evolution of the system for t steps can be obtained from

SSEvolve[rule_, init_, t_, d_Integer] := Nest[FlattenArray[# /. rule, d] &, init, t]

FlattenArray[list_, d_] := Fold[Function[{a, n}, Map[MapThread[Join, #, n] &, a, -{d + 2}]], list, Reverse[Range[d] - 1]]

The analog in 3D of the 2D rule on page 187 is

{1  Array[If[LessEqual[##], 0, 1] &, {2, 2, 2}], 0  Array[0 &, {2, 2, 2}]}

Note that in d dimensions, each black cell must be replaced by at least d + 1 black cells at each step in order to obtain an object that is not restricted to a dimension d - 1 hyperplane.