raytracer.py

quad.fs

#version 130

/* This comes interpolated from the vertex shader */
in vec2 texcoord;

/* The texture we are going to sample */
uniform sampler2D tex;

out vec4 color;

void main(void) {
  /* Well, simply sample the texture */
  color = texture(tex, texcoord);
}

quad.vs

#version 130

/* The position of the vertex as two-dimensional vector */
in vec2 vertex;

/* Write interpolated texture coordinate to fragment shader */
out vec2 texcoord;

void main(void) {
  gl_Position = vec4(vertex, 0.0, 1.0);

  /*
   * Compute texture coordinate by simply
   * interval-mapping from [-1..+1] to [0..1]
   */
  texcoord = vertex * 0.5 + vec2(0.5, 0.5);
}

trace.cs

#version 430 core

layout(binding = 0, rgba32f) uniform image2D framebuffer;

uniform vec3 eye;
uniform vec3 ray00;
uniform vec3 ray01;
uniform vec3 ray10;
uniform vec3 ray11;
uniform vec2 sample_offset;
uniform int frame;

#define MAX_SCENE_BOUNDS    100.0
#define NUM_BOXES           1
#define EPS                 0.000001

const vec3 lightpos = vec3(-1.5, -1.5, 2.5);

const vec3 intensity = vec3(0.8, 0.8, 0.8);
const vec3 ambient = vec3(0.2, 0.2, 0.2);
const vec3 specular = vec3(0.75, 0.75, 0.75);

struct bbox_t {
    float xmin, ymin, zmin;
    float xmax, ymax, zmax;
};

struct mesh_desc_t {
    int vptr;
    int nptr;
    int num_vertices;
    int material_id;
};

struct material_desc_t {
    float Ka, Kd, Ks, Kr, shininess;
};

struct hitinfo_t {
    float t;
    int hit_mesh;
    int hit_vptr;
    int hit_nptr;
};

layout (std430, binding = 1) buffer Materials
{
    material_desc_t materials[];
};

layout (std430, binding = 2) buffer BoundingBoxes
{
    bbox_t box[];
};

layout (std430, binding = 3) buffer MeshDesc
{
    mesh_desc_t desc[];
};

layout (std430, binding = 4) buffer Meshes
{
    float data[];
};

const vec3 diffuse[] = {
    vec3(0.5, 0.5, 0.5),
    vec3(0.2, 0.6, 1.0),
    vec3(1.0, 0.0, 0.0)
};

bool vec3GreaterThan(vec3 a, vec3 b) {
    return (a.x > b.x && a.y > b.y && a.z > b.z);
}

bool intersectBound(vec3 origin, vec3 dir, int i) {
    vec3 vmin = vec3(box[i].xmin, box[i].ymin, box[i].zmin);
    vec3 vmax = vec3(box[i].xmax, box[i].ymax, box[i].zmax);
    if(vec3GreaterThan(origin, vmin) && vec3GreaterThan(vmax, origin))
        return true;
    vec3 tMin = (vmin - origin) / dir;
    vec3 tMax = (vmax - origin) / dir;
    vec3 t1 = min(tMin, tMax);
    vec3 t2 = max(tMin, tMax);
    float tNear = max(max(t1.x, t1.y), t1.z);
    float tFar = min(min(t2.x, t2.y), t2.z);
    return (tNear > 0.0 && tNear < tFar);
}

bool intersectTriangle(vec3 origin, vec3 dir, int ptr, out float dist) {
    vec3 a = vec3(data[ptr], data[ptr + 1], data[ptr + 2]);
    vec3 b = vec3(data[ptr + 3], data[ptr + 4], data[ptr + 5]);
    vec3 c = vec3(data[ptr + 6], data[ptr + 7], data[ptr + 8]);
    vec3 e1 = b - a;
    vec3 e2 = c - a;
    vec3 p = cross(dir, e2);
    float det = dot(e1, p);
    if(abs(det) < EPS) return false;
    det = 1.0 / det;

    vec3 t = origin - a;
    float u = dot(t, p) * det;
    if(u < 0.0 || u > 1.0) return false;
    vec3 q = cross(t, e1);
    float v = dot(dir, q) * det;
    if(v < 0.0 || u + v > 1.0) return false;
    dist = dot(e2, q) * det;
    if(dist > EPS) return true;
    return false;
}

void intersectRay(vec3 O, vec3 D, vec3 P, vec3 Q, out float t, out float s)
{
    t = ((P.x - O.x) * Q.y + (O.y - P.y) * Q.x) / (D.x * Q.y - D.y * Q.x);
    s = (O.x - P.x + t * D.x) / Q.x;
}

vec3 interpolateNormal(vec3 hit_point, hitinfo_t h)
{
    vec3 a = vec3(data[h.hit_vptr], data[h.hit_vptr + 1], data[h.hit_vptr + 2]);
    vec3 b = vec3(data[h.hit_vptr + 3], data[h.hit_vptr + 4], data[h.hit_vptr + 5]);
    vec3 c = vec3(data[h.hit_vptr + 6], data[h.hit_vptr + 7], data[h.hit_vptr + 8]);
    vec3 O = a, D = hit_point - a;
    vec3 P = b, Q = c - b;
    float t, s;
    intersectRay(O, D, P, Q, t, s);

    vec3 Na = vec3(data[h.hit_nptr], data[h.hit_nptr + 1], data[h.hit_nptr + 2]);
    vec3 Nb = vec3(data[h.hit_nptr + 3], data[h.hit_nptr + 4], data[h.hit_nptr + 5]);
    vec3 Nc = vec3(data[h.hit_nptr + 6], data[h.hit_nptr + 7], data[h.hit_nptr + 8]);
    vec3 Nt = normalize(mix(Nb, Nc, s));
    return normalize(mix(Na, Nt, 1 / t));
}

bool isIntersected(vec3 origin, vec3 dir, out hitinfo_t h)
{
    float dist;
    bool hit = false;
    h.t = MAX_SCENE_BOUNDS;
    for(int i = 0; i < desc.length(); i++) {
        if(intersectBound(origin, dir, i)) {
            for(int j = desc[i].vptr; j < desc[i].nptr; j += 9) {
                if(intersectTriangle(origin, dir, j, dist)) {
                    hit = true;
                    if(h.t > dist) {
                        h.t = dist;
                        h.hit_mesh = i;
                        h.hit_vptr = j;
                        h.hit_nptr = j - desc[i].vptr + desc[i].nptr;
                    }
                }
            }
        }
    }
    return hit;
}

const int MAX_TRACE = 4;

struct trace_state_t
{
    vec3 origin;
    vec3 dir;
    vec3 color;
    hitinfo_t h;
    vec3 hit;
    vec3 Ld;
    vec3 N;
};

vec3 trace(vec3 origin, vec3 dir)
{
    trace_state_t stack[MAX_TRACE + 1];
    bool fallback[MAX_TRACE + 1];
    int sp = 0;

    hitinfo_t hl;
    float specular_factor, LdN;

    stack[sp].origin = origin;
    stack[sp].dir = dir;
    stack[sp].color = vec3(0, 0, 0);
    fallback[sp] = false;

    vec3 ambient_contribution;
    vec3 specular_contribution;
    vec3 reflection_contribution;
    vec3 R;

    while(sp >= 0) {
        if(sp >= MAX_TRACE) {
            sp--;
        } else if(fallback[sp]) {
            ambient_contribution = materials[desc[stack[sp].h.hit_mesh].material_id].Ka * ambient;
            reflection_contribution = materials[desc[stack[sp].h.hit_mesh].material_id].Kr * stack[sp + 1].color;
            if(isIntersected(stack[sp].hit, stack[sp].Ld, hl)) {
                stack[sp].color =
                    ambient_contribution + reflection_contribution;
            } else {
                float LdN = dot(stack[sp].Ld, stack[sp].N);
                if(LdN < 0.0) LdN = 0.0;
                R = normalize(2 * LdN * stack[sp].N - stack[sp].Ld);
                specular_factor = dot(R, -stack[sp].dir);
                specular_contribution = vec3(0, 0, 0);
                if(specular_factor > 0) {
                    specular_contribution = materials[desc[stack[sp].h.hit_mesh].material_id].Ks * pow(specular_factor, materials[desc[stack[sp].h.hit_mesh].material_id].shininess) * specular;
                }
                stack[sp].color = ambient_contribution + materials[desc[stack[sp].h.hit_mesh].material_id].Kd * LdN * intensity + specular_contribution + reflection_contribution;
            }
            sp--;
        } else {
            if(isIntersected(stack[sp].origin, stack[sp].dir, stack[sp].h)) {
                stack[sp].hit = stack[sp].origin + stack[sp].h.t * stack[sp].dir;
                stack[sp].Ld = normalize(lightpos - stack[sp].hit);
                stack[sp].N = interpolateNormal(stack[sp].hit, stack[sp].h);
                fallback[sp++] = true;
                stack[sp].origin = stack[sp - 1].hit;
                stack[sp].dir = normalize(stack[sp - 1].dir - 2 * dot(stack[sp - 1].dir, stack[sp - 1].N) * stack[sp - 1].N);
                stack[sp].color = vec3(0, 0, 0);
                fallback[sp] = false;
            } else sp--;
        }

    }

    return stack[0].color;
}

layout (local_size_x = 16, local_size_y = 8) in;
void main(void) {
    ivec2 pix = ivec2(gl_GlobalInvocationID.xy);
    ivec2 size = imageSize(framebuffer);
    if (pix.x >= size.x || pix.y >= size.y) {
        return;
    }
    vec2 pos = vec2(pix.x + sample_offset.x, pix.y + sample_offset.y) / vec2(size.x - 1, size.y - 1);
    vec3 dir = mix(mix(ray00, ray01, pos.y), mix(ray10, ray11, pos.y), pos.x);
    vec3 color = trace(eye, normalize(dir));

    imageStore(framebuffer, pix, mix(imageLoad(framebuffer, pix), vec4(color.rgb, 1.0), 1.0 / float(frame)));
}

tracer.py

from OpenGL.GLUT import *
from OpenGL.GLU import *
from OpenGL.GL import *
from ctypes import c_void_p, c_int, byref
from operator import div, sub, add, mul
from copy import deepcopy

import sys
import numpy as np
import math
import time
import random

NULL = c_void_p(0)

tex = 0
vao = 0
DEFWIDTH, DEFHEIGHT = 600, 720
width, height = 0, 0

# projection / view matrices
viewMatrix = None
projMatrix = None
invMatrix = None

# camera properties
cop = None

# model
meshes = []
materials = []
boundingboxes = []
ssbo_meshes, ssbo_meshesdesc, ssbo_bboxes, ssbo_materials = 0, 0, 0, 0
vssbo, fssbo, nssbo = 0, 0, 0

# shaders
computeProgram = 0
renderProgram = 0

groupsize_x, groupsize_y = 0, 0
sample_lowbound, sample_range = -0.5, 1.0
frame = 0
uni_eye, uni_ray00, uni_ray01, uni_ray10, uni_ray11 = 0, 0, 0, 0, 0
uni_sample_offset, uni_frame = 0, 0

def main():
    glutInit(sys.argv)
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH)
    glutInitWindowSize(DEFWIDTH, DEFHEIGHT)
    glutInitContextVersion(4,2)
    glutInitContextFlags(GLUT_CORE_PROFILE)
    glutCreateWindow('raytracer')
    glutDisplayFunc(display)
    glutReshapeFunc(resize)
    glutIdleFunc(idle)

    resize(DEFWIDTH, DEFHEIGHT)
    readModel()
    addFloor()
    buildSSBOs()
    buildVertexArrays()
    loadShaders()
    initShaders()
    glClearColor(0.,0.,0.,1.)
    #glShadeModel(GL_SMOOTH)
    glEnable(GL_DEPTH_TEST)
    glutMainLoop()
    return

def idle():
    glutPostRedisplay()

def calculateNormalForFace(vertices, normals, f):
    u = np.array(map(sub, vertices[f[1]], vertices[f[0]]))
    v = np.array(map(sub, vertices[f[2]], vertices[f[0]]))
    N = np.cross(u, v).tolist()
    normals[f[0]] = map(add, normals[f[0]], N)
    normals[f[1]] = map(add, normals[f[1]], N)
    normals[f[2]] = map(add, normals[f[2]], N)

def flattenMeshStructure(vertices, faces, normals):
    data_vertices = []
    data_normals = []
    for f in faces:
        for v in f:
            data_vertices.extend(vertices[v])
            data_normals.extend(normals[v])
    return 3 * len(faces), data_vertices + data_normals

def getMeshesData():
    index = 0
    desc, data = [], []
    for i, m in enumerate(meshes):
        desc += [index, 3 * m[0] + index, m[0], m[2]]
        data += m[1]
        index += 6 * m[0]       # 3 for vertex, 3 for normal
    return desc, data

def addFloor():
    global meshes, materials, boundingboxes
    v = [1.0, -1.0, -0.1, 1.0, 1.0, -0.1, -1.0, 1.0, -0.1, 1.0, -1.0, -0.1, -1.0, 1.0, -0.1, -1.0, -1.0, -0.1]
    n = [0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0]
    meshes.append((6, v + n, 1))
    boundingboxes.extend([-1.0, -1.0, -0.1, 1.0, 1.0, 0.1])
    materials += [0.4, 0.0, 0.0, 0.8, 64.0]

def readModel():
    print 'reading model...'
    global meshes, materials, boundingboxes
    vertices, faces = [], []
    min_v = [9999., 9999., 9999.]
    max_v = [-9999., -9999., -9999.]
    with open('body.obj') as f:
        for row in f:
            l = row.split(' ')
            if l[0] == 'v':
                v = map(float, l[1:])
                vertices.append(v)
                min_v = map(min, min_v, v)
                max_v = map(max, max_v, v)
            else:
                faces.append(map(lambda x: int(x) - 1, l[1:]))
    normals = [[0, 0, 0] for i in range(len(vertices))]
    map(lambda f: calculateNormalForFace(vertices, normals, f), faces)
    normals = map(lambda n: map(div, n, [np.linalg.norm(n)] * 3), normals)
    fnormal = open('normal', 'w')
    fnormal.write(repr(normals))
    fnormal.close()

    if min_v[2] < 0:
        delta = 0 - min_v[2]
        vertices = map(lambda v: [v[0], v[1], v[2] + delta], vertices)
        min_v[2] += delta;
        max_v[2] += delta;
    num, data = flattenMeshStructure(vertices, faces, normals)
    materials += [0.8, 0.6, 0.5, 0.05, 64.0]
    meshes.append((num, data, 0))
    boundingboxes.extend(min_v)
    boundingboxes.extend(max_v)

def mapShaderStorageBufferObject(size):
    func = ctypes.pythonapi.PyBuffer_FromMemory
    func.restype = ctypes.py_object
    vp = glMapBuffer(GL_SHADER_STORAGE_BUFFER, GL_READ_ONLY)
    buffer = func(ctypes.c_void_p(vp), size)
    array = np.frombuffer(buffer, dtype=np.dtype(GLfloat))
    return array.tolist()

def createSSBO(array, elemtype, elemsize):
    ssbo = glGenBuffers(1)
    glBindBuffer(GL_SHADER_STORAGE_BUFFER, ssbo)
    arraytype = (elemtype * len(array))
    print 'creating ssbo with actual element list size', len(array)
    glBufferData(GL_SHADER_STORAGE_BUFFER, len(array) * elemsize, arraytype(*array), GL_STATIC_READ);
    glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0)
    return ssbo

def buildSSBOs():
    print 'building SSBO...'
    global ssbo_meshes, ssbo_meshesdesc, ssbo_bboxes, ssbo_materials
    desc, data = getMeshesData()
    ssbo_meshesdesc = createSSBO(desc, GLuint, 4)
    ssbo_meshes = createSSBO(data, GLfloat, 4)
    ssbo_bboxes = createSSBO(boundingboxes, GLfloat, 4)
    ssbo_materials = createSSBO(materials, GLfloat, 4)

def buildTextures():
    print 'building textures...'
    global tex
    if tex:
        glDeleteTextures([tex])
    tex = glGenTextures(1)
    glBindTexture(GL_TEXTURE_2D, tex)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, width, height, 0, GL_RGBA, GL_FLOAT, NULL);
    glBindTexture(GL_TEXTURE_2D, 0)

def buildVertexArrays():
    print 'building vertex arrays...'
    global vao
    vao = glGenVertexArrays(1)
    vbo = glGenBuffers(1)
    glBindVertexArray(vao)
    glBindBuffer(GL_ARRAY_BUFFER, vbo)
    vertices = [-1., -1., 1., -1., 1., 1., 1., 1., -1., 1., -1., -1.]
    arraytype = (GLfloat * len(vertices))
    glBufferData(GL_ARRAY_BUFFER, len(vertices) * 4, arraytype(*vertices), GL_STATIC_DRAW)
    glEnableVertexAttribArray(0)
    glVertexAttribPointer(0, 2, GL_FLOAT, 4, 0, NULL)
    glBindVertexArray(0)

def compileShader(fname, type):
    shader = glCreateShader(type)
    with open(fname, 'r') as f:
        code = f.read()
    glShaderSource(shader, code)
    glCompileShader(shader)
    if glGetShaderiv(shader, GL_COMPILE_STATUS) != GL_TRUE:
        raise RuntimeError(glGetShaderInfoLog(shader))
    return shader

def loadShaders():
    print 'loading shaders...'
    global computeProgram, renderProgram

    quadfs = compileShader('quad.fs', GL_FRAGMENT_SHADER)
    quadvs = compileShader('quad.vs', GL_VERTEX_SHADER)
    tracecs = compileShader('trace.cs', GL_COMPUTE_SHADER)

    computeProgram = glCreateProgram()
    glAttachShader(computeProgram, tracecs)
    glLinkProgram(computeProgram)
    renderProgram = glCreateProgram()
    glAttachShader(renderProgram, quadvs)
    glAttachShader(renderProgram, quadfs)
    glLinkProgram(renderProgram)

def initShaders():
    global groupsize_x, groupsize_y
    global uni_eye, uni_ray00, uni_ray01, uni_ray10, uni_ray11, uni_sample_offset, uni_frame

    glUseProgram(computeProgram)
    groupsize_x, groupsize_y = 16, 8 # glGetProgramiv(computeProgram, GL_COMPUTE_WORK_GROUP_SIZE)
    uni_eye = glGetUniformLocation(computeProgram, 'eye')
    uni_ray00 = glGetUniformLocation(computeProgram, 'ray00')
    uni_ray01 = glGetUniformLocation(computeProgram, 'ray01')
    uni_ray10 = glGetUniformLocation(computeProgram, 'ray10')
    uni_ray11 = glGetUniformLocation(computeProgram, 'ray11')
    uni_sample_offset = glGetUniformLocation(computeProgram, 'sample_offset')
    uni_frame = glGetUniformLocation(computeProgram, 'frame')
    glUseProgram(0)

    glUseProgram(renderProgram)
    tex = glGetUniformLocation(renderProgram, 'tex')
    glUniform1i(tex, 0)
    glUseProgram(0)

def setFrustum(fovy, aspect, near, far):
    global projMatrix
    halfH = math.tan(fovy / 180. * 3.14159 * 0.5) * near
    halfW = aspect * halfH
    projMatrix = np.matrix([
        [near / halfW, 0., 0., 0.],
        [0., near / halfH, 0., 0.],
        [0., 0., (far + near) / (near - far), 2. * far * near / (near - far)],
        [0., 0., -1., 0.]])
    printMatrix('projMatrix', projMatrix)
    #glFrustum(-halfW, halfW, -halfH, halfH, near, far)

def printMatrix(name, m):
    print name,':'
    print '  %.2f\t%.2f\t%.2f\t%.2f\t' % (m[0,0], m[0,1], m[0,2], m[0,3])
    print '  %.2f\t%.2f\t%.2f\t%.2f\t' % (m[1,0], m[1,1], m[1,2], m[1,3])
    print '  %.2f\t%.2f\t%.2f\t%.2f\t' % (m[2,0], m[2,1], m[2,2], m[2,3])
    print '  %.2f\t%.2f\t%.2f\t%.2f\t' % (m[3,0], m[3,1], m[3,2], m[3,3])

def printVector(name, v):
    print name, "%.2f %.2f %.2f" % (v[0], v[1], v[2])

def setCamera(eye, at, up):
    print 'setting camera...'
    global viewMatrix, cop
    cop = deepcopy(eye)
    #gluLookAt(*(eye + at + right))
    eye = np.array(eye)
    at = np.array(at)
    up = np.array(up)
    direction = at - eye
    direction = direction / np.linalg.norm(direction)
    printVector('direction', direction)
    up = up / np.linalg.norm(up)
    right = np.cross(direction, up)
    right = right / np.linalg.norm(right)
    printVector('right', right)
    up = np.cross(right, direction)
    up = up / np.linalg.norm(up)
    printVector('up', up)
    viewMatrix = np.matrix([
        [x for x in right] + [-np.dot(right, cop)],
        [x for x in up] + [-np.dot(up, cop)],
        [-x for x in direction] + [np.dot(direction, cop)],
        [0., 0., 0., 1.0]])
    printMatrix('viewMatrix', viewMatrix)

def calcInvPVMatrix():
    global invMatrix
    invMatrix = np.linalg.inv(projMatrix * viewMatrix)
    printMatrix('inversed PV matrix', invMatrix)

def resize(w, h):
    global width, height
    if w == width and h == height:
        return
    width, height = w, h
    buildTextures()
    glViewport(0, 0, w, h)
    #glMatrixMode(GL_PROJECTION)
    #glLoadIdentity()
    setFrustum(60, float(w) / float(h) if h != 0 else float(w), 1., 2.)
    #glMatrixMode(GL_MODELVIEW)
    eye = [-1.8, -1.8, 0.6]
    at = [0., 0., 0.6]
    right = [0., 0., 1.]
    setCamera(eye, at, right)
    calcInvPVMatrix()

def getEyeRay(x, y, dbg):
    M = np.asarray(invMatrix).tolist()
    x = [x, y, 0., 1.]
    y = []
    if dbg:
        print type(M)
    for row in M:
        if dbg:
            print 'multiplying', row, 'to', x
        y.append(reduce(add, map(mul, row, x)))
    if dbg:
        printVector('y', y)
    y = map(div, y, [y[3]] * 4)
    return map(sub, y[0:3], cop)

def nextPower2(x):
    if x == 0:
        return 1
    x -= 1
    x |= x >> 1
    x |= x >> 2
    x |= x >> 4
    x |= x >> 8
    x |= x >> 16
    x += 1
    return x

def diff(a, b):
    if len(a) != len(b):
        print 'different lengths:', len(a), len(b)
    for i in range(min(len(a), len(b))):
        if abs(a[i] - b[i]) > 0.000001:
            print 'DIFF at position', i

first_time = True
def display():
    global first_time, frame
    frame = frame + 1
    glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
    glUseProgram(computeProgram)

    start_time = time.time()

    # if first_time:
    #     glBindBuffer(GL_SHADER_STORAGE_BUFFER, vssbo)
    #     content = mapShaderStorageBufferObject(len(vertices) * 3 * 4)
    #     flat = []
    #     for v in vertices:
    #         flat.extend(v)
    #     diff(content, flat)
    #     glUnmapBuffer(GL_SHADER_STORAGE_BUFFER)
    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, ssbo_materials)
    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, ssbo_bboxes)
    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, ssbo_meshesdesc)
    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 4, ssbo_meshes)

    glUniform3f(uni_eye, *cop)
    ray = getEyeRay(-1., -1., first_time)
    glUniform3f(uni_ray00, *ray)
    if first_time:
        printVector('ray00', ray)
    ray = getEyeRay(-1., 1., first_time)
    glUniform3f(uni_ray01, *ray)
    if first_time:
        printVector('ray01', ray)
    ray = getEyeRay(1., -1., first_time)
    glUniform3f(uni_ray10, *ray)
    if first_time:
        printVector('ray10', ray)
    ray = getEyeRay(1., 1., first_time)
    glUniform3f(uni_ray11, *ray)
    if first_time:
        printVector('ray11', ray)
    glUniform2f(uni_sample_offset,
        random.random() * sample_range + sample_lowbound,
        random.random() * sample_range + sample_lowbound)
    glUniform1i(uni_frame, frame)
    glBindImageTexture(0, tex, 0, False, 0, GL_WRITE_ONLY, GL_RGBA32F)
    worksize_x, worksize_y = nextPower2(width), nextPower2(height)
    glDispatchCompute(worksize_x / groupsize_x, worksize_y / groupsize_y, 1)

    glBindImageTexture(0, 0, 0, False, 0, GL_READ_WRITE, GL_RGBA32F)
    glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT)
    glUseProgram(renderProgram)
    glBindVertexArray(vao)
    glBindTexture(GL_TEXTURE_2D, tex)
    #if first_time:
    #    with open('d:\\opengl\\texture.dat', 'w') as f:
    #        buf = glGetTexImage(GL_TEXTURE_2D, 0, GL_RGBA, GL_UNSIGNED_BYTE)
    #        f.write(buf)
    #    print len(buf), 'bytes written'
    glDrawArrays(GL_TRIANGLES, 0, 6)
    glBindTexture(GL_TEXTURE_2D, 0)
    glUseProgram(0)
    glutSwapBuffers()
    first_time = False

    duration = time.time() - start_time
    print 'Frame #%d rendering cost %.2f second, estimate fps %.2f' % (frame, duration, 1. / duration if duration > 0.0001 else 0.0)
    return

if __name__ == '__main__':
    main()