Aperiodic Tilings

When you thought that your obsession had run its course and left you in peace, finally, along comes something new and shiny to once more lead you astray. I've written half a dozen posts about sorting networks and genetic algorithms; I thought I had put that topic to rest. What am I to do when I find the Tilings Encyclopedia later the same week? How could I resist this aperiodic fest of self-similar tilings?

The nice thing about the Encyclopedia is not that the tilings are extensively cross-referenced and annotated with the relevant bibliographies; that's all very useful, I don't doubt, for its intended audience; but what I immediately liked was that the tilings are presented as rewriting rules. That should make it easy to replicate them, I thought; I wasn't very wrong, but I find that a functional program is not finished until its structure is strong and regular enough to let it stand on its own legs. After refactoring and restructuring and rewriting the code for a week or so, I had a robust, general program that let me add new tilings, write the tilings in different formats and so on; my last post deals with the details of the back-end mechanism.

Look at the page for the Amman-Beenker rhomb triangle tiling, for instance (go ahead, I'll wait for you). The tiling uses three primitives: two oriented triangles and a parallelogram. Each is "split" or rewritten into a configuration of the same tiles. The process by which the plane is tiled by rewriting is called inflation; the rate of expansion of a primitive tile into a configuration of tiles at each step is called the inflation coefficient.

In typical OCaml fashion, the tyranny of early bound top-level definitions means that I have to sweat the small details first, so that the intended objective is reached in a grand finale by a simple function that synthesizes all the bottom-up development. That's is at odds with writing an engaging literate program, I'm afraid, but it has to be done. The first ingredient, then, in my recipe for pretty pictures is something that lets me put the pieces where they should go. A little bit of linear algebra, in short.

I'll use the most direct representation possible of affine transformations, an affine matrix T such that a vector v gets mapped to v⋅T. Only, I'll be sloppy and identify points on the plane (x, y) ∈ R² with homogeneous vectors (x, y, 1) ∈ A². I always use row vectors and multiply them on the left, which means that the last column of the matrix T is fixed to (0, 0, 1)^T, so I don't need to represent it. To simplify things, I flatten the matrix in column-major order:

type xform = float array

With that, the symmetries of the dihedral group D₂ are simple constant matrices:

let identity = [|  1.;  0.;  0.;  0.;  1.;  0. |]
and fliph    = [|  1.;  0.;  0.;  0.; -1.;  0. |]
and flipv    = [| -1.;  0.;  0.;  0.;  1.;  0. |]
and reflect  = [| -1.;  0.;  0.;  0.; -1.;  0. |]
and qrot     = [|  0.; -1.;  0.;  1.;  0.;  0. |]

The affine transforms are equally direct:

let scale s = [| s; 0.; 0.; 0.; s; 0. |]
and translate (x, y : point) = [| 1.; 0.; x; 0.; 1.; y |]
and skew (x, y : point) = [| 1.; x; 0.; y; 1.; 0. |]
and rotate a =
  let s, c = sin a, cos a in
  [| c; -. s; 0.; s; c; 0. |]

Note that I use homogeneous scaling. The composition of two transforms is simply the product of the corresponding matrices. I denote the composition with the % operator, but perhaps a better choice would have been an infix symbol that highlighted the non-commutative nature of transformations:

let ( % ) (t : xform) (u : xform) : xform =
  [|
    t.(0) *. u.(0) +. t.(3) *. u.(1);
    t.(1) *. u.(0) +. t.(4) *. u.(1);
    t.(2) *. u.(0) +. t.(5) *. u.(1) +. u.(2);
    t.(0) *. u.(3) +. t.(3) *. u.(4);
    t.(1) *. u.(3) +. t.(4) *. u.(4);
    t.(2) *. u.(3) +. t.(5) *. u.(4) +. u.(5);
  |]

To apply a transformation to a point is as simple as vector-matrix multiplication:

let apply (t : xform) (x, y : point) : point =
  (x *. t.(0) +. y *. t.(1) +. t.(2),
   x *. t.(3) +. y *. t.(4) +. t.(5))

The underlying symmetries of the tilings I'll reproduce are polygonal (pentagonal, octagonal, decagonal and so on), so I find easier to work in degrees:

let pi = 3.141592653589793238462643383276
let deg x = x *. 0.017453292519943295769236907684
let rad x = x *. 57.295779513082320876798154814169

The gist of a tiling is captured in the following signature:

module type TILING = sig
  type t
  val name  : string
  val shape : t -> point list
  val color : t -> rgbcolor
  val seed  : (t * xform) list
  val rules : (t * xform) -> (t * xform) list
end

The type t abstracts away the tiles used. name is the full name of the substitution, for purposes of meaningfully labeling the output. shape t returns the polygon corresponding to the tile t as a list of points, and color t returns an RGB triple for that tile. seed is a list of placed tiles with which the inflation process will start. Finally, the meat of the tiling is encoded in rules, that for each tile t returns the expansion as a list of placed tiles. By convention, I'll scale the resulting expansion so that the transformations are local to each tile; in other words, I'll avoid the need to keep track of the overall scale. With that, the core algorithm is captured in the following function:

let iterate f l n =
  let rec go l n =
    if n == 0 then l else
    go (List.concat (List.map f l)) (pred n)
  in go l n

This rather generic function repeatedly expands a tiling l by applying n times the rules f. Interpreted under the list monad, iterate is nothing more than functional iteration of the monadic function f; but clearly bringing the entire machinery of the non-deterministic monad for just this is overkill.

With this, I have all the necessary pieces to actually implement the drawing algorithm. In keeping with the theme of the previous post, I'll output EPS code to a temporary file which I'll then open with Preview. The Generator is parameterized by a TILING:

module Generator (S : TILING) = struct

But, there's a snag to solve: as I've set up things, the final size of the tiling will be that of the seed, and I use, in general, shapes with sides of length 1. This will make the final tiling the size of a pixel! Precisely this is the drive behind using "abstract" back-ends: I can draw the tiling once to get its "real" size, rescale it, and then draw it for real:

  module Tiling (G : GRAP) = struct
    module M = Monad(G.M)
    open M.Op

    let draw = M.mapm_ (fun (t, x) ->
      let vs = List.map (apply x) (S.shape t) in
         G.color (S.color t)
      >> G.poly ~fill:true vs
      >> G.gray 0.
      >> G.poly ~fill:false vs)
  end

Remember that a GRAP back-end contains the monad M that implements the back-end. Drawing a tiling involves taking each placed tile (t, x), transforming the polygon S.shape t by the placement transform x, and drawing it. Again, this is completely generic on the "drawing" method. I take advantage of this genericity first to find the bounding box of the tiling:

  let bounds tiling =
    let module T = Tiling(Bounds) in
    Bounds.bounds (T.draw tiling)

Then, I draw the tiling simply by setting an appropriate line width and using the correct back-end:

  let show width height tiling =
    let module T = Tiling(Eps) in
    Eps.show width height (T.M.seq (Eps.weight 0.5) (T.draw tiling))

This pair of functions is the pay-off of so much refactoring and abstraction! Finally, generate ties all together by finding the bounding box, scaling the tiling and outputting it:

  let transform x' = List.map (fun (t, x) -> (t, x % x'))

  let generate dim levels =
    let tiling = iterate S.rules S.seed levels in
    let bbox   = bounds tiling in
    let width  = bbox.Bounds.right -. bbox.Bounds.left
    and height = bbox.Bounds.bottom -. bbox.Bounds.top in
    let size   = dim /. width in
    let xform  = translate (-. bbox.Bounds.left, -. bbox.Bounds.top) % scale size in
    show (width *. size) (height *. size) (transform xform tiling)
end

This function has a rather arbitrary design decision: whether to scale to the bounding box of the seed, or to that of the final tiling. If the expansion were confined to the boundary of the parent tile, it would make no difference (except for the work involved!); not all tilings, however, are so tidy. One of the most famous aperiodic tilings, Penrose's, overstep the boundaries of the parent tile.

The job is done, except for the actual tilings. First, the Amman-Beenker rhomb-triangle tiling:

module AmmanBeenkerRhombTriangle : TILING = struct
  type t = Tu | Td | Pa

  let name = "Amman-Beenker"

  let rho  = sqrt 2.
  and rho' = sqrt 0.5

  let shape = function
  | Tu -> [0., 0.; rho', rho'; rho, 0.; 0., 0.]
  | Td -> [0., 0.; rho, 0.; rho', -. rho'; 0., 0.]
  | Pa -> [0., 0.; 1., 0.; 1. +. rho', -. rho'; rho', -. rho'; 0., 0.]

  let color = function
  | Tu -> (0.976, 0.4, 0.063)
  | Td -> (0.973, 0.796, 0.4)
  | Pa -> (0.804, 0.796, 0.855)

  let seed = 
    let xform = rotate (deg 45.) in
    [Tu, xform; Td, xform]

  let rules (t, x) =
    let x' = scale (rho -. 1.) % x in match t with
  | Tu -> [
    Tu, rotate (deg (-135.)) % translate (1., 1.) % x';
    Pa, rotate (deg (-90.)) % translate (1. +. rho', 1. +. rho') % x';
    Tu, rotate (deg 135.) % translate (2. +. rho', rho') % x';
    Pa, rotate (deg 0.) % translate (1. +. rho', rho') % x';
    Td, rotate (deg 180.) % translate (1. +. rho, 0.) % x'; 
  ]
  | Td -> [
    Td, rotate (deg 135.) % translate (1., -1.) % x';
    Pa, rotate (deg 135.) % translate (1. +. rho', -1. -. rho') % x';
    Td, rotate (deg (-135.)) % translate (2. +. rho', -. rho') % x';
    Pa, rotate (deg 45.) % translate (1. +. rho', -. rho') % x';
    Tu, rotate (deg 180.) % translate (1. +. rho, 0.) % x'; 
  ]
  | Pa -> [
    Pa, x';
    Td, rotate (deg 0.) % translate (1., 0.) % x';
    Tu, rotate (deg 135.) % translate (2. +. rho, -1.) % x';
    Pa, rotate (deg 0.) % translate (1. +. rho, -1.) % x';
    Td, rotate (deg 180.) % translate (1. +. rho +. rho', -1. -. rho') % x';
    Tu, rotate (deg (-45.)) % translate (rho', -. rho') % x';
    Pa, rotate (deg (-90.)) % translate (1. +. rho, 0.) % x';
  ]
end

There's not much to say about this code, except that it encodes precisely the data presented in its Encyclopedia page. Let me just add that I found encoding the transformations quite error-prone; but then again I'm not a visual thinker. The result is:

let () =
  let module G = Generator(AmmanBeenkerRhombTriangle) in
  G.generate 256. 4

This tile could make excellent backgrounds; note that the weaving of triangles and parallelograms creates a somewhat disturbing optical illusion. Another beautiful tiling is the Robinson triangle tiling:

module RobinsonTriangle : TILING = struct
  type t = Su | Sd | Au | Ad

  let name = "Robinson.Triangle"

  let phi  = 0.25 *. (sqrt 5. +. 1.)
  let phi' = 0.25 *. (sqrt 5. -. 1.)
  and eta  = sqrt (0.125 *. (5. +. sqrt 5.))
  and eta' = sqrt (0.125 *. (5. -. sqrt 5.))

  let shape = function
  | Su -> [0., 0.; 2. *. phi , 0.; phi ,    eta'; 0., 0.]
  | Sd -> [0., 0.; 2. *. phi , 0.; phi , -. eta'; 0., 0.]
  | Au -> [0., 0.; 2. *. phi', 0.; phi',    eta ; 0., 0.]
  | Ad -> [0., 0.; 2. *. phi', 0.; phi', -. eta ; 0., 0.]

  let color = function
  | Su -> (0., 0., 0.376)
  | Sd -> (0.808, 0.796, 0.89)
  | Au -> (1., 0.8, 0.4)
  | Ad -> (1., 0.4, 0.)

  let seed =
    let t = rotate (deg 36.) in
    [ Su, t; Sd, t ]

  let rules (t, x) =
    let x' = scale (0.5 /. phi) % x in match t with
  | Su -> [
    Su, rotate (deg (-144.)) % translate (0.5 +. phi, eta) % x';
    Ad, rotate (deg ( -36.)) % translate (0.5 +. phi, eta) % x';
    Sd, reflect % translate (1. +. 2. *. phi, 0.) % x';
  ]
  | Sd -> [
    Sd, rotate (deg 144.) % translate (0.5 +. phi, -. eta) % x';
    Au, rotate (deg  36.) % translate (0.5 +. phi, -. eta) % x';
    Su, reflect % translate (1. +. 2. *. phi, 0.) % x';
  ]
  | Au -> [
    Au, rotate (deg 108.) % translate (1., 0.) % x';
    Su, rotate (deg (-108.)) % translate (0.5, 2. *. eta *. phi) % x';
  ]
  | Ad -> [
    Ad, rotate (deg (-108.)) % translate (1., 0.) % x';
    Sd, rotate (deg 108.) % translate (0.5, -2. *. eta *. phi) % x';
  ]
end

The seed brings together to obtuse triangles to form a parallelogram; there is, however, no way to "cut" the tile into a rectangular tile with matching borders, so it won't work as a background:

let () =
  let module G = Generator(RobinsonTriangle) in
  G.generate 256. 6

Finally, the classic Penrose rhomb tiling:

module PenroseRhomb : TILING = struct
  type t = R0 | R1

  let name = "Penrose.Rhomb"

  let phi  = 0.25 *. (sqrt 5. +. 1.)
  let phi' = 0.25 *. (sqrt 5. -. 1.)
  and eta  = sqrt (0.125 *. (5. +. sqrt 5.))
  and eta' = sqrt (0.125 *. (5. -. sqrt 5.))

  let shape = function
  | R0 -> [0., 0.; eta', -. phi ; 0., -2. *. phi ; -. eta', -. phi ; 0., 0.]
  | R1 -> [0., 0.; eta , -. phi'; 0., -2. *. phi'; -. eta , -. phi'; 0., 0.]

  let color = function
  | R0 -> (0.973, 0.796, 0.4)
  | R1 -> (0.08, 0.027, 0.4)

  let seed = [
    R0, identity;
    R0, rotate (deg 72.);
    R0, rotate (deg 144.);
    R0, rotate (deg 216.);
    R0, rotate (deg 288.);
  ]

  let rules (t, x) =
    let x' = scale (2. *. phi') % x in match t with
  | R0 -> [
    R0, reflect % translate (0., -2. *. phi) % x';
    R1, rotate (deg 36.) % translate (eta', -. phi) % x';
    R0, rotate (deg 144.) % translate (0., -1. -. 2. *. phi) % x';
    R0, rotate (deg (-144.)) % translate (0., -1. -. 2. *. phi) % x';
    R1, rotate (deg (-36.)) % translate (-. eta', -. phi) % x';
  ]
  | R1 -> [
    R1, rotate (deg 108.) % translate (-. eta', phi' -. 0.5) % x';
    R1, rotate (deg (-108.)) % translate (eta', phi' -. 0.5) % x';
    R0, rotate (deg 108.) % translate (0., -1.) % x';
    R0, rotate (deg (-108.)) % translate (0., -1.) % x';
  ]
end

As you can see, the tiles overflow their parents. I used a "sun" seed to generate a pentagonally-symmetrical result:

let () =
  let module G = Generator(PenroseRhomb) in
  G.generate 256. 4

I've uploaded the complete program, with the code in the last post and in this, to the OCamlCore Forge. There are a number of really nice tilings I'll maybe some day implement noted in comments at the end of the file. Or, find the ones you like and add a snippet to the forge. Have fun!