Art with Algorithms¶
On some approaches and algorithms that are usefule in procedural art.
- Grammars
- Chaotic Maps
- Noise
- Putting it all together
Grammars - Context-Free¶
L-Systems¶
L-Systems are named for Aristid Lindenmayer, the biologist who first introduced them in 1968. He devised them as a tool for modelling plant growth. They are loosely equivalent to context-free grammars, with a few differences. Most notably, in a L-System production rules are expanded in parallel per evaluation rather than in-series.
$G_L = (\Sigma, S, \delta)$
$\Sigma$ - alphabet
$S$ - start symbol
$\delta$ - production rules
Example: Lindemayer's Algae¶
In [11]:
def l_system(state, rules):
'''Generate strings for an L-System.
state -- a string with the system's initial state
rules -- a dictionary of production rules
'''
while True:
next_state = ""
for sym in state:
# check if production should be applied
if sym in rules:
next_state += rules[sym]
# else copy verbatim
else:
next_state += sym
# save state and yield
state = next_state
yield state
# a simple system of Lindenmayer's
algae = l_system('A', {'A':'AB', 'B':'A'})
# print the first seven strings
for state in zip(range(7), algae):
print(*state)
Example: Fractal Tree¶
In [13]:
tree = l_system('0', {'1':'11', '0': '1[-0]+0'})
for i, state in zip(range(8), tree):
if i < 4: print(i, state)
reset(ctx)
ctx.set_line_width(2)
ctx.translate(256, 512)
ctx.rotate(-pi/2)
draw_state(ctx, next(tree), angle=pi/4)
ctx.stroke()
In [14]:
display.Image(srf.write_to_png())
Out[14]:
Example: Sierpinski Triangle¶
In [15]:
# the sierpinski triangle
tri = l_system('A+B+B', {'A':'A+B-A-B+A', 'B':'BB'})
for _ in range(8):
next(tri)
# setup to draw
reset(ctx)
ctx.translate(484, 512)
ctx.rotate(-pi/2)
ctx.set_line_width(0.5)
# draw and display
draw_state(ctx, next(tri), angle=2*pi/3)
ctx.stroke()
In [16]:
display.Image(srf.write_to_png())
Out[16]:
Random L-System¶
In [193]:
def rand_str(alpha, low, high):
'''Return a random string'''
from random import choice, randint
return ''.join(choice(alpha) for _ in range(randint(low, high)))
# the alphabet
alpha = 'ABC-+'
# the initial state
state = rand_str('ACB', 1, 5)
# the production rules
rules = {x: rand_str(alpha, 1, 8) for x in 'ABC'}
# a random system
rand_sys = l_system(state, rules)
print('state: {}\nrules: {}'.format(state, rules))
In [201]:
display.Image(srf.write_to_png())
Out[201]:
In [214]:
display.Image(srf.write_to_png())
Out[214]:
In [215]:
def l_system_csg(state, rules):
while True:
next_state = ""
while state:
# check if any production should be applied
skip = False
for key in rules:
match = state.find(key)
if match == 0:
skip = True
next_state += rules[key]
state = str(state[len(key):])
break
# else copy verbatim
if not skip:
next_state += state[0]
state = str(state[1:])
# save state and yield
state = next_state
yield state
In [216]:
thing_gen = l_system_csg('A', {'+B':'BB-A', 'B':'A', 'A':'A+BB'})
for i, state in zip(range(9), thing_gen):
if i < 6: print(i, state)
reset(ctx)
ctx.identity_matrix()
ctx.translate(256, 256)
ctx.rotate(-pi/2)
ctx.set_line_width(0.5)
draw_state(ctx, next(thing_gen), angle=pi/4)
ctx.stroke()
In [217]:
display.Image(srf.write_to_png())
Out[217]:
Grammars - Recursively Enumerable¶
Turmites¶
A Turmite is a Turing Machine run on a 2D tape that posesses the additional internal state of an "orientation". Of course like an n-tape or n-head Turing Machine, Turmites are strictly equivalent in capability to their single-tape cousings. Turmites work quite simply:
- Read a symbol on the tape.
- Write a symbol onto the tape.
- Move in one of four directions.
In [218]:
import numpy as np
def turmite(tape, symbols=10, states=5):
from numpy.random import randint
state = 0 # internal state
bearing = 0 # direction turmite is facing
ind = tuple(randint(0, x, 1)[0] for x in tape.shape)
# Tape Symbols States
# ------------ --------------------------------------
# 0 (write sym, move %3, next state) (...)
# 1 (write sym, move %3, next state) (...)
# 2 (write sym, move %3, next state) (...)
dims = (symbols, states)
# randomize the delta function
delta = np.zeros((*dims, 3), dtype=np.uint8)
delta[:, :, 0] = randint(0, symbols, dims) # symbol to write
delta[:, :, 1] = randint(0, 4, dims) # move direction
delta[:, :, 2] = randint(0, states, dims) # next state
while True:
# read tape
write, move, state = delta[tape[ind], state]
bearing = (bearing + move) % 4
# write tape symbol
tape[ind] = write
# move tape head
offset = ((0, 1), (1, 0), (-1, 0), (0, -1))[bearing]
ind = tuple(np.mod(np.add(ind, offset), tape.shape))
yield
In [257]:
anim
Out[257]:
In [241]:
# setup for drawing
fig = pylab.figure()
fig.set_figwidth(15)
# setup a blank tape
tape = np.zeros((100, 100), dtype=np.uint8)
# run six turmites
for n in range(6):
ax = fig.add_subplot(2, 3, n+1)
ax.set_axis_off()
np.ndarray.fill(tape, 0)
tm1 = turmite(tape, states=n+2)
for _ in range(10000): next(tm1)
pyplot.imshow(tape, vmin=0, vmax=9)
Chaotic Maps¶
Chaotic maps are discrete-time dynamic systems. Many are derived from hamiltonian systems or other physical systems. The following Ikeda map is an example Map, used to model dielectric insulators.
In [234]:
# Ikeda Map
def ikeda(ctx, scale=512, u=0.9, n=100, steps=100):
from math import cos, sin
# select 'n' random points
points = np.random.random((n, 3))
points[:, 0] *= scale
points[:, 1] *= scale
for _ in range(steps):
for i in range(n):
x, y, t = points[i]
ctx.move_to(x, y)
x_new = 1.0 + u*(x*cos(t) - y*sin(t))
y_new = u*(x*sin(t) + y*cos(t))
t_new = 0.4 - 6.0/(1.0 + x**2 + y**2)
ctx.line_to(x_new, y_new)
ctx.stroke()
points[i] = x_new, y_new, t_new
# prepare canvas
ctx.set_line_width(0.1)
ctx.set_source_rgba(1, 1, 1, 0.05)
reset(ctx, color=(0, 0, 0, 1))
ctx.translate(256, 256)
ctx.scale(10)
ikeda(ctx, n=1000)
In [235]:
display.Image(srf.write_to_png())
Out[235]:
In [57]:
reset(ctx, color=(0, 0, 0, 1))
ctx.translate(256, 256)
ctx.scale(10)
ikeda(ctx, u=0.75, n=2000)
display.Image(srf.write_to_png())
Out[57]:
In [244]:
def ds(rp=2, base=.5, wrap=(False, False)):
import numpy
power = 1
data = None
while True:
# length of side
side = 2**power + 1
sz = [side, side]
# create new numpy image
image = numpy.zeros(sz)
# if wrapped on either (x,y)
if wrap[0]:
sz[1] -= 1
if wrap[1]:
sz[0] -= 1
# copy & expand last image
if type(data) != type(None):
for x in range(data.shape[0]):
for y in range(data.shape[1]):
image[(2*x) % sz[0], (2*y) % sz[0]] = data[x, y]
# set corners to random data
else:
for x in ((0, 0), (sz[0]-1, 0), (0, sz[1]-1), (sz[0]-1, sz[1]-1)):
image[x] = max(min(image[x] + (numpy.random.random() - .5) + base, 1.0), 0.0)
# generate new pixels
for x in range(1, side, 2):
for y in range(1, side, 2):
xp1 = (x + 1) % sz[0]
yp1 = (y + 1) % sz[1]
mod = (numpy.random.random() - 0.5)/rp**power
image[x-1, y] = mod + (image[x-1, yp1] + image[x-1, y-1]) / 2
image[xp1, y] = mod + (image[xp1, yp1] + image[xp1, y-1]) / 2
image[x, y-1] = mod + (image[x-1, y-1] + image[xp1, y-1]) / 2
image[x, yp1] = mod + (image[x-1, yp1] + image[xp1, yp1]) / 2
image[x, y] = mod + (image[x-1, y-1] + image[x-1, yp1] + image[xp1, y-1] + image[xp1, yp1]) / 4
# normalize all pixels
for x in numpy.nditer(image, op_flags=['readwrite']):
x[...] = x if x >= 0. else 0.
x[...] = x if x <= 1. else 1.
# save current image
data = image[0:sz[0], 0:sz[1]]
power += 1
# expand to size (side, side)
if wrap[1]:
image[sz[0]] = image[0]
if wrap[0]:
image[..., sz[1]] = image[..., 0]
# return image
yield image
In [248]:
anim
Out[248]:
Example: Diamond-Square¶
In [71]:
# diamond square
def dsn(n, *args):
p_gen = ds(*args)
for it, s in zip(p_gen, range(n)):
square = it
return p_gen, square
p_gen, planet = dsn(7, 1.8, .5, (True, True))
pyplot.imshow(cm.Vega20c(planet))
Out[71]:
In [75]:
# show potential ds output
fig = pylab.figure()
fig.set_figwidth(15)
for n in range(0,6):
ax = fig.add_subplot(2, 3, n+1)
ax.set_axis_off()
tmp_gen = ds(1.8, .5, (False, False))
for it, s in zip(tmp_gen, range(6)):
copy = it
pyplot.imshow(cm.Vega20c(copy))
In [76]:
def shadow(image, a=.5, b=0, frac=.05, tol=.25, below=True):
b *= image.shape[0]
tol *= image.shape[0]
for x in range(image.shape[0]):
for y in range(image.shape[1]):
fx = a*x + b
if y - fx < tol and y > fx:
image[x, y, 0:3] *= (y - fx)/tol * (1 - frac) + frac
elif (below and y <= fx) or (not below and y >= fx):
image[x, y, 0:3] *= frac
In [77]:
copy = cm.Vega20c(planet.copy())
shadow(copy)
pyplot.imshow(copy)
Out[77]:
In [78]:
def lens_distort(image, alpha = []):
from math import pi as π
from math import atan2, cos, sin, asin, hypot
# copy image
copy = image.copy()
# side
side = image.shape[0]
for cx in range(side):
for cy in range(side):
# normalized coordinates
nx, ny = cx/side - 0.5, cy/side - 0.5
# normalized visual distance
nr = hypot(nx, ny)
if nr > 0.5:
image[cx, cy] = [0., 0., 0., 0.]
continue
# angle of x,y
Ï• = atan2(ny, nx)
# distorted distance
d = asin(2*nr) / π
# map pixel on new image -> pixel on old image
image[cx, cy] = copy[round(side*(0.5 + d * cos(Ï•))), round(side*(0.5 + d * sin(Ï•)))]
In [80]:
copy = cm.Vega20c(planet.copy())
grid(copy)
fig = pylab.figure()
fig.set_figwidth(15)
fig.add_subplot(1, 2, 1)
pyplot.imshow(copy)
fig.add_subplot(1, 2, 2)
lens_distort(copy)
pyplot.imshow(copy)
pylab.show()
In [88]:
copy = cm.Vega20c(planet.copy())
gifify((roll(copy),), len(copy))
Out[88]:
In [97]:
c_gen = ds(1.5, .25, (True, True))
for it, s in zip(c_gen, range(7)):
clouds = it
clouds = cm.gray(clouds)
pyplot.imshow(clouds)
for x in range(clouds.shape[0]):
for y in range(clouds.shape[1]):
clouds[x, y] = [1., 1., 1., clouds[x, y, 0]]
In [99]:
copy = cm.Vega20c(planet.copy())
gen = (roll(copy, axis=(0,1)), roll(clouds, -2, 0), roll(clouds))
gifify(gen, 2*len(clouds), shadow, lens_distort)
Out[99]: