graphics.py

import numpy

DIMENSIONS = 3
SCREEN_DIMENSIONS = 2

class Palette(object):
    def __init__(self, zero, one):
        self.zero, self.one = numpy.array(zero), numpy.array(one)

    def color(self, c):
        assert 0. <= c <= 1.
        precise = (1. - c) * self.one + c * self.zero
        integer = numpy.rint((precise * 255)).astype(int)
        return tuple(integer)
PALETTE_BLUE = Palette([0., 0., 1.], [1., 1., 1.])
assert PALETTE_BLUE.color(0) == tuple([255, 255, 255])
assert PALETTE_BLUE.color(1) == tuple([0, 0, 255])
assert PALETTE_BLUE.color(.5) == (128, 128, 255)


def create_identity():
    return numpy.identity(DIMENSIONS + 1)

def create_projection():
    """
    Returns a 4x4 projection matrix that, applied to 3d homogenous
    coordinates [x, y, z, 1] of a point in world-space, would give
    [x'w, y'w, z'w, w] such that x' = x'w/w and y' = y'w/w are the
    point's screen-space coordinates.
    
    @see project()
    """
    # TODO
    pass
assert numpy.allclose(
    create_projection(),
    create_projection().dot(create_projection()))

def create_scaling(axes):
    """
    When called with len(axes) = N, returns a matrix that scales each
    axis in N-dimensional homogenous coordiates by the given amount.
    
    e.g., to zoom in the Y axis in 3D by 2, use create_scale([1, 2, 1]).
    """
    # TODO
    pass
assert numpy.allclose(
    create_scaling([1, 2, 3]).dot([1, 1, 1, 1])[:-1],
    [1, 2, 3])

def create_translation(point):
    """
    When called with len(point) = N, returns a matrix that translates
    N-dimensional homogenous coordiates by the given amount.
    
    e.g., to move in the Y axis in 2D by 2, use create_translation([0, 2]).
    """
    # TODO
    pass
assert numpy.allclose(
    create_translation([1, 2, 3]).dot([1, 1, 1, 1])[:-1],
    [2, 3, 4])

def create_rotation(angle, axis1, axis2):
    """
    Returns a DIMENSIONS-dimensional rotation matrix that rotates points
    by angle on the axis1-axis2 plane.
    
    e.g., to rotate by 90 degrees in the XZ plane ([1, 0, 0] -> [0, 0, 1]),
    use create_rotation(pi/2, 0, 2).
    """
    # TODO
    pass
assert numpy.allclose(
    create_rotation(numpy.pi/2, 0, 1).dot([1, 0, 0, 1])[:-1],
    [0, 1, 0])
assert numpy.allclose(
    create_rotation(numpy.pi/2, 1, 2).dot([1, 0, 0, 1])[:-1],
    [1, 0, 0])
assert numpy.allclose(
    create_rotation(numpy.pi/4, 0, 1).dot([1, 0, 0, 1])[:-1],
    [1/numpy.sqrt(2), 1/numpy.sqrt(2), 0])

def create_cube(dimensions):
    """
    Returns a `dimensions`-dimensional cube, e.g. a square (2D) or a hypercube (4D).
    The cube spans from [0, ..., 0] to [1, ..., 1].
    
    The shape is returned as a tuple (vertices, lines) where vertices is a KxD table
    (each column is a point in D dimensions, and there are K points) and lines is a
    list of pairs of indices of rows in vertices.
    
    e.g., a 1D cube (a line) would be
    ([[0],
      [1]],
     [(0, 1)]),
    and a 2D cube (a square) would be
    ([[0, 1, 0, 1],
      [0, 0, 1, 1]],
     [(0, 1), (1, 2), (2, 3), (3, 0)]),
    which describes the points: {A: (0, 0), B: (1, 0), C: (0, 1), D: (1, 1)}
    and line segments between the points (A, B), (B, C), (C, D) and (D, A),
    or alternatively the "less orderly"
    ([[0, 1, 0, 1],
      [0, 1, 1, 0]],
     [(2, 1), (3, 0), (0, 2), (1, 3)]),
    which describes the points: {A: (0, 0), B: (1, 1), C: (0, 1), D: (1, 0)}
    and line segments between the points (C, B), (D, A), (A, C) and (B, D).
    """
    # TODO
assert numpy.allclose(
    create_cube(2)[0],
    [[0, 1, 0, 1],
     [0, 0, 1, 1]])
assert numpy.allclose(
    create_cube(2)[1],
    [(0, 1), (2, 3), (0, 2), (1, 3)])
assert create_cube(3)[0].shape[1] == 8
assert create_cube(3)[1].shape[0] == 12

def project(projection, camera, source):
    """Apply camera matrix, then projection matrix, to source, which
    has 3D points as columns.
    """
    # add a row of 1 values to the bottom, to convert 3D coordinates
    # to 3D homogenous coordinates
    homogenous = numpy.vstack((source,
                               numpy.ones((1, source.shape[1]),
                                          dtype=source.dtype)))
    # apply camera, then projection
    nonhomogenous = (projection.dot(camera.dot(homogenous)))
    # translate homogenous to 2D points: divide x, y by w
    vertices = nonhomogenous[:-1, :]
    vertices[:2] /= nonhomogenous[-1, :]
    assert vertices.shape == (SCREEN_DIMENSIONS + 1, source.shape[1])
    return vertices

render.py

import pyglet
from pyglet.gl import *

import graphics
from graphics import create_projection, create_translation, create_scaling, create_identity,\
    create_rotation, create_cube, project, PALETTE_BLUE

DIMENSIONS = 3

LINE_WIDTH = 2.0

CAPTION = "4D"
X_SIZE, Y_SIZE = 640, 480
INTERVAL = 1.0 / 60.0


class Renderer(object):
    def __init__(self):
        self.x_size = X_SIZE
        self.y_size = Y_SIZE
        self.window = pyglet.window.Window(self.x_size, self.y_size,
                                           caption=CAPTION, vsync=False)
        self.projection = create_projection()
        self.projection = (
            create_translation((300., 300., 0.))).dot(
                self.projection).dot(
                    create_scaling((100., 100., 1.)))
        self.camera = create_identity()

        self.create_bounding_cube()

        pyglet.clock.schedule_interval(self.update, INTERVAL)
        self.fps_counter = pyglet.clock.ClockDisplay()
        self.window.event(self.on_draw)

    def create_bounding_cube(self):
        bounding_vertices, bounding_lines = create_cube(DIMENSIONS)
        bounding_vertices *= .96
        bounding_vertices += .2
        self.bounding_cube = (bounding_vertices, bounding_lines)

    def on_draw(self):
        "Automaticaly called upon every clock tick, renders the board"
        self.window.clear()
        glEnable(GL_BLEND)
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
	glMatrixMode(GL_MODELVIEW)
	glLoadIdentity()

        glDisable(GL_TEXTURE_2D)
        self.render(*self.bounding_cube)
        glEnable(GL_TEXTURE_2D)

        self.fps_counter.draw()

    CAMERA_ROTATION = (
        create_translation([.5, .5, .5]).dot(
        create_rotation(.003, 0, 1)).dot(
        create_translation([-.5, -.5, -.5]))
    )
    def update(self, dt):
        self.camera = self.camera.dot(self.CAMERA_ROTATION)

    def render(self, vertices, lines, palette=PALETTE_BLUE):
        glLineWidth(LINE_WIDTH)
        vertices = project(self.projection, self.camera, vertices)
        for a, b in lines:
            x1, y1, c1 = vertices[:, a]
            x2, y2, c2 = vertices[:, b]
            if c1 < 0 or c2 < 0:
                continue
            c1 = min(1., c1)
            c2 = min(1., c2)
            pyglet.graphics.draw(2, GL_LINES,
                                 ('v2f', (x1, y1, x2, y2)),
                                 ('c3B', palette.color(c1) + palette.color(c2)))


if __name__ == '__main__':
    renderer = Renderer()
    pyglet.app.run()