Algorithms for rendering quasicrystal tilings

Quasicrystals are aperiodic spatial tilings, and 2D quasicrystals have been used for centuries in Islamic art. I find these patterns beautiful, and I think more people should know how to render them.

You can render 2D quasicrystals as a sum of plane waves. They can also be constructed using aperiodic tilings (e.g. 4-fold Ammann-Beenker tiles and 5-fold Penrose Tiles), and there are fractal constructions that recursively subdivide a set of tiles using a smaller set of the same tiles (e.g. 5- and 7-fold tilings). Many more constructions can be found on the Tilings Encyclopedia

This post focuses on the "cut and project" construction, which is easy to implement in code. The routines shown here can be found in this iPython notebook on Github.

The cut-and-project construction

The simplest way to render crisp, geometric quasicrystals is the "cut and project" construction. This views the quasicrystal as a two-dimensional projection of a planar slice through a N-dimensional hyper-lattice. Let's implement this algorithm. 

A (hyper) lattice

First, we need to define the idea of a hyperlattice. This is a N-dimensional generalization of a cubic crystal lattice. It tiles ND space using ND hypercubes. We can represent the center of each hypercube (hypercell? hyperpixel?) as a N-D vector of integers
$$
\mathbf q_0 = [ x_1, ... x_N ],\,\,x_i\in\mathbb Z
$$
The associated N-cube, then, includes all points on a unit hypercube $\left[-\tfrac 1 2, \tfrac 1 2\right]^N$ offset to location $\mathbf q_0$:
$$
\mathbf q = \left[-\tfrac 1 2, \tfrac 1 2\right]^N + \mathbf q_0
$$

An irrational projection onto a 2D plane

To construct the quasicrystal, we cut through this ND lattice with a 2D plane. To get a aperiodic slice, this plane must not align with any periodic direction in the hyperlattice. 

A good choice is to project each direction of the hyperlattice onto an evenly-spaced set of basis directions, centered at the origin. We define a 2D $(s,t)$ plane $\mathcal P$ in terms of two N-D basis vectors, $\mathbf u, \mathbf v \in \mathbb R^N$.

$$
\begin{aligned}
\mathbf p &= \mathbf u \cdot s + \mathbf v \cdot t,\,\,s,t\in\mathbb R
\\
\mathbf u &= [ \cos(\theta_0), \dots , \cos(\theta_{N-1}) ]
\\
\mathbf v &= [ \sin(\theta_0), \dots , \sin(\theta_{N-1}) ]
\\
\theta_i &= \frac {\pi\cdot i}{N}
\end{aligned}
$$

To render the quasicrystal, cut the N-D lattice along this plane, and project the exposed surface of this cut isometrically onto the 2D (s,t) plane $\mathcal P$. I've searched for an elegant algorithm for this to no avail. It seems like there should be some generalization of e.g. Bresenham's line algorithm to finding the hypercubes (hyperpixels?) associated with the cut. But, inelegant solutions work: it suffices to find the set of N-cubes in the ND lattice that intersect this plane.

$$
\mathcal Q = \left\{ \mathbf q_0 \bigm| \mathcal P \cap \mathbf q\ne \varnothing \right\}
$$

 

 

Algorithm 1: brute-force (wrong, but easy)

An easy to code (but terrible) way to find points in $\mathcal Q$ is by a brute-force search: Check lots of points in the s,t plane, and quantize them to the nearest vertex $\mathbf q_0$ of the hyper-lattice:
$$
\mathcal Q = \left\{ \mathbf q_0 \bigm| \exists (s,t)\in\mathbb R^2 \text{ s.t. } \left \lfloor \mathbf u s + \mathbf v t \right \rceil =  \mathbf q_0 \right\},
$$
where $\lfloor\cdot\rceil$ denotes rounding to the nearest integer vector. This is massively inefficient: many $(s,t)$ points map to the same lattice point $\mathbf q$, and it will miss points if the grid resolution is too coarse. But, for small renderings and run-once code, it does the job. Here it is Python3 code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

from scipy.spatial.distance import squareform, pdist from pylab import * from numpy import * N = 4 # Dimension of crystal θ = linspace(0,pi,N+1)[1:] # Arrange basis vectors evenly on [0,π) A = array([cos(θ),sin(θ)]) # Convert from complex notation L = 501 # Number of grid points to search over dx = 5 # Size of square portion of plane to check # Generate the search grid x = linspace(-dx,dx,L) xy = array(meshgrid(x,x)).reshape(2,L*L) # Collapse all grid points to the nearest hyper-cell p = unique(np.round(A.T@xy),axis=1) def plotpoints(p): u = pinv(A).T@p D = squareform(pdist(p.T)) e = {tuple(sorted(e)) for e in zip(*where(abs(D-1)<1e-6))} xy = concatenate([[u[:,a],u[:,b],[NaN,NaN]] for a,b in e]) plot(*xy.T,lw=0.5,color='k') axis("equal") axis("off") plotpoints(p)

Figure 1: rendered output of a rhobic tiling for the N=4 quasicrystal, computed by brute-force search of the cut-and-project approach.

Testing for intersections

A better way is to start with a point that you know touches the plane. Then, "grow" the crystal from this point by checking if adjacent sites on the lattice also intersect the cutting plane. The plane defined above always intersects the N-cube at $\mathbf q=[0,..,0]$, so this is a good place to start.

This mathematics stack exchange question gives two ways to check for cube-plane intersection. They include (a) using linear programming, and (b) checking directly based on a geometric solution. Both are slower than I'd like, and I suspect there is an even faster/more elegant solution out there. 

Algorithm 2: test with linear programming

As per this answer, we can check for plane-cube intersection using Linear Programming (LP). Mirko Stojiljković's "Hands-On Linear Programming: Optimization With Python" is a good, accessible introduction to linear programming in Python. 

The code below searches for N-cubes intersecting the cutting plane iteratively. It includes optimizations, like memoizing the intersection tests, and exploiting the N-fold symmetry to recycle previous test results.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

from scipy.optimize import linprog # Number of search iterations Ndepth = 10 # Directions in which to search D = int32(concatenate([eye(N),-eye(N)])) # Matrix to rotate a lattice point π/N in 2D plane S = zeros((N,N)) S[ -1,0 ] = -1 S[:-1,1:] = eye(N-1) # Matrix encoding LHS of linear inequalities for LP Aub = concatenate([ A.T,-A.T]) # Use linear programming to see if plane # intersects hypercube ±0.5 of p def check_lp(q): return linprog(c=[1,1],A_ub=Aub, b_ub=concatenate([.5+q,.5-q]), method="revised simplex", options={"tol":1e-9}).success # Cache results to save time (memoization) cache = {} def checkcrystal_lp(q): # Convert test point to immutable tuple for cache key k = tuple(int32(q)) if not k in cache: # Recompute intersection test if not in cache # Use symmetry: reduce all tests to points in 1st sector h = angle([1,1j]@A@q) cache[k]=check_lp(q) if 0<=h<=(pi/N*1.1) else checkcrystal_lp(S@q) return cache[k] # Start with a seed point at zero and build outwards Q,allQ = {(0,)*N},set() # Iteratively search in each direction to add points for i in range(Ndepth): Q = {tuple(q) for p in Q for q in D+p if checkcrystal_lp(q)}-allQ allQ |= Q figure(figsize=(10,10)) plotpoints(array(list(allQ)).T)

Figure 2: Rhombic quasicrystal tiling for N=4, computed via iterative search. Linear programming (algorithm 2) was used in this case, but other tests for plane-hypercube intersection work similarly (see  algorithms 3 and 4).

Algorithm 3: testing using a geometric solution

This answer outlines another way to test for plane-cube intersection. It involves projecting the problem onto a space where the 2D plane collapses to a point at the origin, and checking whether the shadow of the N-cube in this subspace contains this origin point. 

There are a few ways to do this. For example, you could calculate the convex hull of the N-2 projection, and test if it contains the origin. The code below uses another approach: It checks whether any of the N-2 sub-facets of the projected N-cube contains the origin, stopping early if it finds one.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68

from scipy.linalg import null_space from itertools import product, combinations # Null space: collapses cutting plane to point at origin C = null_space(A).T # All vertices of an N-cube O = array([*product(*([0,1],)*N)])-0.5 # All vertices of a 2-cube (used to generate N-2 sub-facets) F = array([*product(*([0,1],)*2)]) # A single corner ("origin" vertex) of an N-2 cube # plus N-2 vectors defining the other points on the N-2 cube # Used to test if sub-facet contains the origin i = [0]+[*2**arange(N-2)] G = array([*product(*([0,1],)*(N-2))])[i,:] # All possible pairs of dimension for generating # all possible N-2 sub-facets of the N-cube I = [*combinations(range(N),2)] # Given a list of N-cube vertices, construct a list # of all N-2-cube sub-facets. subfacets = [] for ij in I: # For every pair of directions, and for every 2-cube # in each pair of directions... Generate all possible # N-2-cube sub-facets. v0 = F@eye(N)[ij,:] v1 = G@eye(N)[[*{*range(N)}-{*ij}]] subfacets.extend(v0[:,None,:]+v1[None,:,:]) vxid = int32(subfacets)@2**arange(N,dtype='i') # Check if a sub-facet contains the origin def subfacet_contains(v): p0 = v[0] vk = v[1:,:]-p0 q = -np.linalg.solve(vk@vk.T,vk@p0) return all(q>=0) and all(q<=1) # Check if N-2 projection of N-cube contains the origin # (thereby checking for plane-cube intersection) def in_plane(q): # Generate list of all N-2 dimensional sub-facets # of unit N-cube centered at q p = (q+O).T sf = (C@p)[:,vxid].transpose(1,2,0) for v in sf: if subfacet_contains(v): return True return False cache = {} def checkcrystal2(p): k = tuple(int32(p)) if not k in cache: h = angle([1,1j]@A@p) s = in_plane(p) if 0<=h<=(pi/N*1.1) else checkcrystal2(S@p) cache[k]=s return cache[k] Q,allQ = {(0,)*N},set() for i in range(Ndepth): Q = {tuple(q) for p in Q for q in D+p if checkcrystal2(q)}-allQ allQ |= Q figure(figsize=(10,10)) plotpoints(array(list(allQ)).T)

Algorithm 4: Intersection of $n$ half-planes in $\mathcal O [ n\log(n)]$

This algorithm is harder to code, but gives fast performance that scales favorably as the dimension of the crystal grows. To understand the geometry, let's return to the brute-force algorithm. 

In the image on the below, we've chosen a random color for each cell in an $N=5$ hyperlattice. We've iterated over the points in the cutting plane, and colored them according to which cell they reside in. This is easy to compute: we simply use the defined basis $A$ to project each point on the plane into $N$ dimensions, then round the coordinates to the nearest integer

The image on the right was rendered by taking $N$ planar square-waves, and summing them together modulo 2. This illustrates the close connection between the quasicrystal as a sum of $N$ plane waves, and the quasicrystal as a diagonal slices through an $N$-dimensional hyperlattice.

Theses plane waves cut the plane into strips. Regions where $N$ strips overlap in the plane correspond to cells in the hyperlattice that intersect the cutting plane. Thus, we can reduce the test for plane-hypercube intersection to a test for whether the $N$ strips, oriented in different directions, and offset based on the position of the hypercube, have a nonempty intersection. 

To do this, first, transform the problem into a set of $2N$ linear inequality constraints. Each constraint defines a half-plane. The plane intersects the hypercube if the intersection of these half-planes is nonempty. So far, this is identical to the linear programming approach of Algorithm 2.

There are fast algorithms for testing whether the intersection of $n$ half-planes is nonempty. This approach by Preparata and Muller does it in $n\log(n)$. These notes by Dave Mount explain the geometry underlying these algorithms. 

The approach is: 

• First, identify any vertical bounding lines. These demarcate spans of the $x$ axis of the plane. Their intersection $x\in[x_0,x_1]$ defines the bounds for a search procedure (below). (We can always choose an orientation such that one of the crystal axes is vertical.)
• Split the remaining lines into those that bound the half-plane from below, and those that bound the half plane from above, where "below" and "above" are defined in terms of the $y$ coordinate of the plane.
• These two groups of lines define and upper and lower feasibility regions, which are convex. The boundary of these regions can be interpreted as curves $y_l(x)$ and $y_u(x)$
• Test if the the intersection between the upper/lower feasibility regions is nonempty by finding the minimum of $y(x) = y_l(x)-y_u(x)$ on $x\in[x_0,x_1]$
• If any $\exists x\in[x_0,x_1]\text{ s.t. } y(x)<0 $, then the plane intersects the hypercube.
• This can be checked via binary search, looking for the point where $t(s)$ changes sign, and stopping early if any point satisfying all the constraints is found.

Figure 4: Illustration of $\mathcal O(n \log(n))$ test for plane-hypercube intersection. The intersection problem is reduces to testing whether the intersection of $N$ strips in the plane is nonempty. The problem is simplified by letting one of these strips be vertical, restricting the search region to a portion of the $x$ axis. Each remaining strip is expressed as the intersection between an upper-bounding and lower-bounding half-plane. The intersection of all upper-bounding half-planes (black) forms an upper bounding curve (blue). Similarly for the lower bound (orange). If there is any $x$ such that lower-upper is negative, then the intersection is nonempty. This can be checked via binary search.

This intersection test can be used with the iterative search routines shown for Algorithms 2 and 3. There is a demonstration in this iPython notebook. The Python3 source code for the binary search is below.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

def check_pm(p,eps=1e-3,maxiter=1000): p = int32(p) # Get slopes and spacing of boundary lines # This assumes you've set things up so that θ[-1]=π m = 1/tan(θ[:-1]) d = 1/sin(θ[:-1]) # Evaluate upper/lower bounds lp,up = p[:-1]-0.5, p[:-1]+0.5 def feasibility_surface(p,x): u,l = x*m+up*d, x*m+lp*d y = np.max(l)-np.min(u) dy = m[np.argmax(l)]-m[np.argmin(u)] return y,dy # Initial bounds and feasability checks # If either bound satisfies, success # If better point not in bounds, fail x0,x1 = p[-1]-0.5, p[-1]+0.5 y0,dy0 = feasibility_surface(p,x0) y1,dy1 = feasibility_surface(p,x1) if y0 < 0 or y1 < 0: return True if dy0>=0 or dy1<=0: return False for i in range(maxiter): # Guess min by following slope at bounds xm = (y1-y0-dy1*x1+dy0*x0)/(dy0-dy1) assert xm>=x0-eps assert xm<=x1+eps # Midpoint works? still feasible? ym,dym = feasibility_surface(p,xm) if ym<0: return True if dym==0: return False # Binary search: move inward if not done if dym<0: # If better point not in bounds, fail if dym>=0 or abs(x0-xm)<eps: return False x0,y0,dy0 = xm,ym,dym else: if dym<=0 or abs(x1-xm)<eps: return False x1,y1,dy1 = xm,ym,dym # Should have returned succes/fail before getting here raise RuntimeError("Maximum iterations exceeded: p=%s"%p)

---

  Create 

These approaches generate (a subset of) the vertices of a quasicrystaline tiling. The routines shown here can be found in this iPython notebook on Github (edit: the $\mathcal O(n\log(n))$ algorithm is here). You can connect these points to render a rhombic tiling, or calculate mosaic-like Voronoi cells. For the most part, I just render the tilings and then color them manually. The tiles can be subdivided, or used as a structure to support detailed elaborations. Make something beautiful (: