Random noise textures on the GPU in WebGL

Procedural rendering routines often need pseudorandom numbers. For graphics rendering, we usually prefer correct-ish and fast code over properly "correct" code. My approach to random number generation is usually to intermittently seed a weak pseudo-random number generator, like a linear-feedback shift register, from a stronger source like Javascript's Math.random(). 

On the GPU, things are not so simple: we want a pseudorandom algorithm that is fast, local, and parallel. For compatibility, we should also have this RNG use an 8-bit RGBA texture for its internal state. Bitops are a bit tricky in WebGL shaders (though not impossible), so we'll use a linear congruential generator.

Pixel-art-style Conway's game-of-life using a tile shader in WebGL

Implementation of Conway's game of life cellular automata with a retro-style shader. Added as example on webgpgpu Github.

[view in browser]

[smaller version for mobile]

Sketch of WebGL implementation:

• Underlying cellular automata runs Conway's game of life
• States as game tiles: new cells: small blue; live cells: blue; died: skulls
• Life diffuses out coloring nearby areas
• A single RGBA texture stores the game state. We render to/from this texture as a framebuffer.
• Each color channel in [0,1] can be mapped to a [0,255] byte value
• The green channel stores the game of life state
• The blue channel stores the diffusing field
• The red channel outputs a tile ID value that is passed to a tile shader to render the game

Retro-styled forest-fire cellular automata in WebGL

Added as example on webgpgpu Github.

This is a cellular automata analogue of the model that we used for developmental retinal waves in this paper. 

[view in browser]

Key rules

• Vegetation "grows" occasionally (flip a coin to increment amount of trees)
• Fires have a small chance of starting spontaneously
• Intensity of ignited fires decay down to zero, at which point they leave ash
• Fires spread to adjacent cells with a probability that depends on their intensity
• Probability of a cell igniting given a fire is related to how much vegetation there is

Sketch of WebGL implementation:

• A single RGBA texture stores the game state
• Each color channel in [0,1] can be mapped to a [0,255] byte value
• A "noise" texture creates pseudorandom numbers
• The green channel stores the level of vegetation (BURNT,NONE,GRASS,SCRUB,TREE)
• The blue channel stores the intensity of the fire
• The red channel outputs a tile ID value that is passed to a tile shader to render the game

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.

Box Drawing Characters + Marching Squares = Marching Boxes?

I couldn't figure out the name of this to search for it online, and it wasn't mentioned on Wikipedia's marching squares page, so… here's to reinventing wheels ( :

Sebastian Blomfield's "Lord of the Manor" game inspired me to explore retro-game engine coding. So far, I've got a basic tile-shader that renders tiled 2D game environments (e.g. game of life, forest fire simulation). Next up: what if I want to draw maps where the adjacent tiles of e.g. walls are continuous?

Box-drawing characters are a nice way to generate rectilinear paths in a tile-based game environment. By checking whether the tiles above/below/left/right are also walls, we can select one of 16 box drawing symbols to create continuous walls:

Note that this extends the usual unicode box-drawing characters with five additional tiles: one for an isolated (disconnected) wall, and four for "dead ends". These suffice for rendering square hallways, e.g:

With some tweaks to the tiles, we can coax box-drawing tiles to generate more curved and diagonal paths:

See here for an example page that renders dendriform growth using the above tile sets (refresh to get a new randomly-selected tile set).

Marching squares is another nice algorithm for rounding off corners on a map (e.g.). This doesn't quite work out-of-the-box for our tile shader, however, since marching squares accepts the values for the corners of each tile, and I want it to operate on the tile values directly. It also uses 16 tiles:

So! box-drawing is good for lines, and marching squares is good for borders/contours. Can we combine these and get a best-of-both-worlds approach? One that renders tiles based on their current value and the values of their neighbors?

Yes. I'm not sure what its name is, so I'm calling it "marching boxes" for now. Marching-boxes extends the box-drawing tiles with 31 tiles (below). These provide versions of the corner, T-junction, and cross-junction tiles that account for the possibility of nearby filled-in regions. The resulting tiles can render contours around filled-in regions as well as sharp lines.

Including background and foreground tiles, marching boxes uses 48 tiles.

This could be nice for any pattern that has a mixture of lines and enclosed regions, like a dungeon connected by hallways, or a water system that includes both lakes and rivers.

Click here to view a WebGL demo of marching-boxes tile rendering

Visualizing z-transforms in the the complex plane

These notes (pdf) , along with the script plotting_helper.py , are available as an ipython notebook on Github .

Polar and rectangular coordinates

Recall the rectangular $z=x+ iy$ and polar $z=r e^{i\theta}$ forms of a complex number.

Coloring the points in the complex plane

The two maps below are colored representations of the magnitude (left) and phase (right) in polar coordinates. Observe that the phase rotates $2\pi$ around zero. To connect this to the analysis of poles and zeros in our z-domain transfer function, we can think of these maps as coloring the identity map. For example, it could be a z-domain transfer function of the form: $$ U(z) = z $$ This map has has one zero (at zero), and one pole (at infinity).

The way that phase behaves near 0 and $\infty$ in the identity map resembles the behavior of phase around zeros and poles in more complicated s-domain and z-domain transfer functions. In the identity map, going around the origin counterclockwise adds $2\pi$ to the phase. This is equivalent to encircling infinity clockwise . In a transfer function, each zero (pole) can be viewed as a shifted (for poles: and inverted) version of the identity map.

This allows us to determine the # of poles - # zeros by taking a closed contour integral in the complex plane, which can help us assess the stability of a sytem. For example, in continuous time we need there to be no poles in the right half plane. In discrete time, we need to check that there are no poles outside the unit circle.

This follows from Cauchy's argument principle and is important intution for understanding the Nyquist stability criterion in both the continuous (s-domain) and discrete (z-domain). However, we do not need to go through the proof of the argument principle to understand it intuitively.

Maps in the complex plane

By coloring points in the complex plane, we can visualize functions from $\mathbb C\to\mathbb C$

A shift by 1 (zero at -1)

Note that the $2\pi$ phase rotation around the zero, located at -1

Squaring (two zeros at zero)

Note that the phase rotation doubles, and is now $4\pi$ rather than $2\pi$ around zero.

Inverting (pole at 0)

Note that the phase rotation has reversed.

Let's make an unstable conjugate pair z-domain transfer function

$$ U(z) = \frac 1 {(1 + \alpha e^{i\pi/4} z^{-1})(1 + \alpha e^{-i\pi/4} z^{-1})},\,\, \alpha < 1 $$

Phase rotation vanishes around regions that contain equal # zeros and poles

Poles and zeros cancelling inside unit circle: Bode and Nyquist plots

Phase and magnitude together

We can plot the output of a function as a shaded color, where the hue represents the phase and the brightness the absolute magnitude. (We'll use this in a bit to visualize transormation in the complex plan in a single plot)

Just for fun: Fractals as iterated polynomial maps

Just for fun: the iterated map $z_n = z_n^2 + c$ can generate Julia set fractals. Here's an example for $c = -0.81 + 0.05i$, iterated 7 times.

Approximate conversions from Laplace to Z domain

Forward

$$ z = e^{sT} \approx 1 + sT \Rightarrow sT \approx z-1 $$

Backward

$$ 1/z = e^{-sT} \approx 1 - sT \Rightarrow sT \approx \frac{z-1}{z} $$

Tustin/Bilinear

$$ z = \frac {e^{\tfrac 1 2 sT}}{e^{-\tfrac 1 2 sT}} \approx \frac {1 + \tfrac 1 2 sT} {1 - \tfrac 1 2 sT} \Rightarrow sT \approx 2\frac {z-1}{z + 1} $$

Applying each of these transformations in the complex plane...

Applied to a Laplace domain transfer function w. two conjugate pair poles

• All transfomations introduce some error
• Forward Euler can be unstable
• Backward Euler
  • Maps the point at $\infty$ to 0
  • Preserves stability
  • Can distort frequences and phase response
• Tustin transformation
  • Maps the point at $\infty$ to -1
  • Preserves stability
  • Can distort frequences and phase responses
• Transforming poles and zeros
  • Poorly approximates the linear time-invariant frequency response

Comparing frequency domain responses

Note that, even through the system attained via forward Euler is unstable, it still appears to have a reasonable frequency response. This is because we cannot tell whether the poles are inside or outside the unit circle from the frequency response alone, and highlights the importance of checking the stability before interpreting freqeuency reponse plots.

( In these plots, we use an arbitrary (normalized) frequency, and set the sampling rate of the discrete system to be twice this unit frequency (T=0.5) )

Extra example 1

$$ \frac{10^{-4} z^3}{(z-1)(1.15z^2-2.05z+1)} $$

Extra example 2

$$ \frac 1 {z^2(z-1)} $$