Generates a Sierpinski triangle

flame.py

from PIL import Image
import math
from array import array
from collections import namedtuple
from functools import partial
import random
import operator

I = namedtuple('I', ['w', 'h', 'd'])

def i_mk(wh, t='i', d=3): # image make
    w, h = wh
    im = I(w, h, [array(t)]*d)
    for k in range(w*h):
        for l in range(d):
            im.d[l].append(0)
    return im

def i_g(img, rc): # image get
    r, c = rc
    return tuple(map(lambda img_ch: img_ch[img.w*r + c], img.d))

def i_s(img, rc, color): # image set
    r, c = rc
    for i, img_channel in enumerate(img.d):
        img_channel[img.w*r + c] = color[i]

def draw(pil_img, img): # draw an I to a PIL.Image
    r_max, c_max = pil_img.size
    for r in range(r_max):
         for c in range(c_max):
             v = i_g(img, (r,c))
             pil_img.putpixel( (r,c), v)
    return pil_img

def rnd_weight(choices): # choose at random weighted from [(item, prob), ...]
    space = {}
    current = 0
    for choice, weight in choices:
        if weight > 0:
            space[current] = choice
            current += weight
    rand = random.uniform(0, current)
    for key in sorted(space.keys() + [current]):
        if rand < key:
            return choice
        choice = space[key]
    return None

def f0(x): return (x[0]/2.0,
                   x[1]/2.0)
def f1(x): return ((x[0]+.1)/2.0,
                   x[1]/2.0)
def f2(x): return (x[0]/2.0,
                   (x[1]+.1)/2.0)

funs = [
    (f0, 1),
    (f1, 1),
    (f2, 1)
]

def ifsi(funs, x): # ifs iteration
    f = rnd_weight(funs)
    res = f(x)
    return res

def fp(f, x0, N): # fixed point iterate f up to N iterations
    x = x0
    for i in range(N):
        x = f(x)
    return x

def view_t(x, s, t):
    return map(lambda x, t: s*x+t, x, t)

def quant(x, sz):
    return map(lambda x, s: int(s*(x/2.0+0.5)),
               x, sz)

def IFS(funs, x0, N):
    return fp(partial(ifsi, funs), x0, N)

def tonemap(v, m):
    if m == 0: return (0, 0, 0)
    v = v/float(m)
    v = int(255*math.log(1+v)/math.log(2))
    return (v, )*3

def main():
    size = (256, 256)
    pil_img = Image.new('RGB', size)
    img = i_mk(size)
    h = i_mk(size, d=1)
    N = 10000
    M = 20

    h_m = 0
    for i in range(N): # build a histogram
        x = (random.uniform(-1, 1), random.uniform(-1, 1))
        x_ = IFS(funs, x, M)
        x_ = view_t(x_, 13, (-0.5, -0.5))

        x_ = quant(x_, (img.w, img.h))
        hx_ = i_g(h, x_)
        if hx_[0] > h_m: h_m = hx_[0]
        i_s(h, x_, (hx_[0]+1,)) # incr. histogram

    for r in range(size[0]): # map histogram to image
        for c in range(size[1]):
            t = tonemap(i_g(h, (r, c))[0], h_m)
            i_s(img, (r, c), t)

    pil_img = draw(pil_img, img)
    pil_img.save("foo.png", "PNG")

if __name__ == '__main__':
    main()