inspired by http://www.reddit.com/r/mathpics/comments/ookpa/i_decided_to_play_with_the_chaos_game_after_that/ spins, grows and then changes over time

chaoscube.py

#!/usr/bin/python

# inspired by
# http://www.reddit.com/r/mathpics/comments/ookpa/i_decided_to_play_with_the_chaos_game_after_that/

# also, not actually a cube

from OpenGL.GL import *
from OpenGL.arrays import vbo
from OpenGL.GL.shaders import *
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *

import numpy
import numpy.linalg as linalg
import random
import math
from threading import Thread

import Image

class ImageWriter(Thread):
    def __init__(self, image, filename):
        Thread.__init__(self)
        self.image = image
        self.filename = filename
    def run(self):
        self.image.save(self.filename)

aspect = 1.0

tetrahedron_corners = numpy.array([ [1,1,1],
                                    [-1,-1,1],
                                    [-1,1,-1],
                                    [1,-1,-1] ], dtype=numpy.float32)

cube_corners = numpy.array([ [ 1., 1., 1.],
                             [ 1., 1.,-1.],
                             [ 1.,-1., 1.],
                             [ 1.,-1.,-1.],
                             [-1., 1., 1.],
                             [-1., 1.,-1.],
                             [-1.,-1., 1.],
                             [-1.,-1.,-1.] ], dtype=numpy.float32)

twocube_corners = numpy.array([ [ 1., 1., 1.],
                                [ 1., 1.,-1.],
                                [ 1.,-1., 1.],
                                [ 1.,-1.,-1.],
                                [-1., 1., 1.],
                                [-1., 1.,-1.],
                                [-1.,-1., 1.],
                                [-1.,-1.,-1.],
                                [ 2., 2., 2.],
                                [ 2., 2.,-2.],
                                [ 2.,-2., 2.],
                                [ 2.,-2.,-2.],
                                [-2., 2., 2.],
                                [-2., 2.,-2.],
                                [-2.,-2., 2.],
                                [-2.,-2.,-2.],
                                ], dtype=numpy.float32)

phi = (1.0+math.sqrt(5))/2.0
raw_icosahedron_corners = []
for (a,b) in [(1,1), (1,-1), (-1,1), (-1,-1)]:
    raw_icosahedron_corners.append( [0.0, a*1, b*phi] )
    raw_icosahedron_corners.append( [a*1, b*phi, 0] )
    raw_icosahedron_corners.append( [a*phi, 0, b*1] )
icosahedron_corners = numpy.array(raw_icosahedron_corners, dtype=numpy.float32)

theta = 0.0
frame = 0
def Render():
    global theta
    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0, aspect, 0.5, 1000.0)

    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()

    #theta = (glutGet(GLUT_ELAPSED_TIME) / 60.0) % 360.0
    #theta /= (180.0 / math.pi)
    theta += math.pi / 180.0

    gluLookAt(20.0 * math.sin(theta), 20.0 * math.cos(theta), 20.0,
              0, 0, 0,
              0, 0, 1)

    #glEnable(GL_NORMALIZE)
    glEnable(GL_BLEND)

    glClearColor(0,0,0.1,1)
    glClear(GL_COLOR_BUFFER_BIT)

    if True:
        #glUseProgram(shader)
        vbobj.bind()
        glEnableClientState(GL_VERTEX_ARRAY)
        glVertexPointer(3, GL_FLOAT, 0, vbobj )

        cbobj.bind()
        glEnableClientState(GL_COLOR_ARRAY)
        glColorPointer(4, GL_FLOAT, 0, cbobj )

        glDrawArrays(GL_POINTS, 0, vdata.shape[0])

        cbobj.unbind()
        glDisableClientState(GL_COLOR_ARRAY)

        vbobj.unbind()
        glDisableClientState(GL_VERTEX_ARRAY)
        #glUseProgram(0)

    glutSwapBuffers()

    if os.access("frames", os.F_OK):
        buffer=glReadPixels(0,0,1920,1080,GL_RGB,GL_UNSIGNED_BYTE)
        im=Image.fromstring("RGB", (1920,1080), buffer)
        global frame
        ImageWriter(im, "frames/frame-%06i.jpg" % (frame,)).start()
        frame += 1

blocksize = 512
vblock = numpy.zeros((blocksize,3), dtype=numpy.float32)
cblock = numpy.zeros((blocksize,4), dtype=numpy.float32)

def idle(*args):
    global point
    global idx

    scaling = random.uniform(0.5, 1.5)
    for x in range(0, blocksize):
        corner = corners[ random.randint(0, corners.shape[0]-1) ]
        direction = corner - point
        direction *= scaling
        point += direction
        vblock[x][0:3] = point
        cblock[x][3] = 0.1
        foo = direction / linalg.norm(direction)
        foo = (foo - numpy.min(foo))
        foo /= numpy.max(foo)
        cblock[x][0:3] = foo / linalg.norm(foo)
        #block[x][3] = scaling - 0.5
    idx += blocksize
    start = idx % vdata.shape[0]
    vbobj[slice(start, start+blocksize)] = vblock
    cbobj[slice(start, start+blocksize)] = cblock

    glutPostRedisplay()

def update(*args):
    glutTimerFunc(40, update, 0)
    #glutPostRedisplay()

def key(*args):
    if args[0] == '\x1b':
        glutDestroyWindow(window_id)
        try:
            sys.exit(0)
            out = open("foo.matrix","w")
            for pt in range(0, vdata.shape[0]):
                out.write("%f %f %f\n" % tuple(vdata[pt]))
            out.close()
        finally:
            sys.exit(0)

def reshape(width, height):
    global aspect
    aspect = float(width)/float(height) if (height>0) else 1.0
    glViewport(0, 0, width, height)

if __name__ == "__main__":
    glutInit([])

    glutInitDisplayString("rgba>=8 depth>16 double")
    glutInitWindowSize(1920, 1080)
    window_id = glutCreateWindow("chaos cube")

    glutDisplayFunc(Render)
    glutReshapeFunc(reshape)
    glutKeyboardFunc(key)
    glutIdleFunc(idle)

    #glUseProgram(0)

    glutTimerFunc(40, update, 0)

    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
    glEnable(GL_MULTISAMPLE)

    corners = tetrahedron_corners
    corners *= 3.5

    vdata = numpy.zeros((2**20, 3))
    vdata = numpy.array(vdata, dtype=numpy.float32)
    vdata = 2.0 * (vdata - 0.5)

    cdata = numpy.zeros((2**20, 4))
    cdata = numpy.array(cdata, dtype=numpy.float32)

    point = numpy.array([ random.random(), random.random(), random.random() ],
                        dtype=numpy.float32)

    idx = 0

    VERTEX_SHADER = compileShader("""
void main() {
  gl_Position = gl_ModelViewProjectionMatrix * vec4(gl_Vertex.xyz, 1.0);
}
""", GL_VERTEX_SHADER)
    
    FRAGMENT_SHADER = compileShader("""
void main() {
 gl_FragColor = vec4(1.0,0.7,0.7,0.1);
}
""", GL_FRAGMENT_SHADER)

    shader = compileProgram(VERTEX_SHADER, FRAGMENT_SHADER)

    vbobj = vbo.VBO(vdata, usage=GL_STREAM_DRAW)
    cbobj = vbo.VBO(cdata, usage=GL_STREAM_DRAW)

    glutMainLoop()