raytrace.c

// clang -g -I ../IntenseLogic/src raytrace.c -o raytrace -lGL -lGLEW -lSDL2 -L ../IntenseLogic/build -lilcommon -lilmath -lilutil -lilgraphics -ldl -lm -pthread

#include <GL/glew.h>
#include <SDL2/SDL.h>
#include <stdio.h>
#include <stdbool.h>
#include <sys/time.h>
#include <math.h>
#include <assert.h>
#include <pthread.h>

#include "graphics/glutil.h"
#include "math/vector.h"
#include "math/quaternion.h"
#include "common/world.h"
#include "util/array.h"
#include "util/loader.h"

const unsigned int max_depth = 3;

il_world *world;

struct material {
    il_vec3 colour;
    il_vec3 diffuse;
    il_vec3 specular;
    float specular_co;
};

struct sphere {
    il_positionable pos;
    struct material mat;
    float radius;
};
IL_ARRAY(struct sphere,) spheres;

struct box {
    il_positionable pos;
    struct material mat;
};
IL_ARRAY(struct box,) boxes;

struct light {
    il_positionable pos;
    il_vec3 colour;
    float radius;
};
IL_ARRAY(struct light,) lights;

struct ray {
    unsigned int id, root;
    unsigned short depth, x, y;
    il_vec3 origin;
    il_vec3 dir;
    il_vec3 opacity;
};
struct ray_chunk {
    struct ray_chunk *next;
    unsigned length;
    struct ray rays[10000];
} *rays_head, *rays_tail, *rays_free;
unsigned ray_counter = 1;
unsigned num_rays_free;
pthread_mutex_t ray_mutex;

struct ray_chunk *ray_chunk_new()
{
    struct ray_chunk *chunk;
    if (rays_free) {
        chunk = rays_free;
        rays_free = rays_free->next; // pop off the head of the free list
        --num_rays_free;
    } else {
        chunk = malloc(sizeof(struct ray_chunk));
    }
    chunk->next = NULL;
    chunk->length = 0;
    return chunk;
}

void ray_chunk_free(struct ray_chunk *chunk)
{
    if (num_rays_free < 100) {
        free(chunk); // don't let the free list get too big
        return;
    }
    chunk->next = rays_free;
    rays_free = chunk;
    ++num_rays_free;
}

struct ray *emit_ray(struct ray r)
{
    pthread_mutex_lock(&ray_mutex);
    if (rays_tail->length >= 10000) {
        struct ray_chunk *chunk = ray_chunk_new();
        rays_tail->next = chunk;
        rays_tail = chunk;
    }
    size_t id = rays_tail->length++;
    rays_tail->rays[id] = r;
    pthread_mutex_unlock(&ray_mutex);
    return &rays_tail->rays[id];
}

struct diffuse {
    unsigned root;
    unsigned short depth, x, y;
    il_vec3 origin;
    il_vec3 opacity;
};
struct diffuse_chunk {
    struct diffuse_chunk *next;
    size_t length;
    struct diffuse diffuse[10000];
} *diffuse_head, *diffuse_tail, *diffuse_free;
unsigned int num_diffuse_free;
pthread_mutex_t diffuse_mutex;

struct diffuse_chunk *diffuse_chunk_new()
{
    struct diffuse_chunk *chunk;
    if (diffuse_free) {
        chunk = diffuse_free;
        diffuse_free = diffuse_free->next; // pop off the head of the free list
        --num_diffuse_free;
    } else {
        chunk = malloc(sizeof(struct diffuse_chunk));
    }
    chunk->next = NULL;
    chunk->length = 0;
    return chunk;
}

void diffuse_chunk_free(struct diffuse_chunk *chunk)
{
    if (num_diffuse_free < 100) {
        free(chunk); // don't let the free list get too big
        return;
    }
    chunk->next = diffuse_free;
    diffuse_free = chunk;
    ++num_diffuse_free;
}

struct diffuse *emit_diffuse(struct diffuse r)
{
    pthread_mutex_lock(&diffuse_mutex);
    if (diffuse_tail->length >= 10000) {
        struct diffuse_chunk *chunk = diffuse_chunk_new();
        diffuse_tail->next = chunk;
        diffuse_tail = chunk;
    }
    size_t id = diffuse_tail->length++;
    diffuse_tail->diffuse[id] = r;
    pthread_mutex_unlock(&diffuse_mutex);
    return &diffuse_tail->diffuse[id];
}

struct link {
    unsigned int id, root;
    unsigned short depth, x, y;
    il_vec3 opacity;
    il_vec3 colour;
};
IL_ARRAY(struct link,) links;
struct link root_links[800*600];

struct contact {
    struct material *mat;
    float dist;
    il_vec3 normal;
};

bool test_sphere(struct sphere *sphere, struct ray *ray, struct contact *out)
{
    // http://wiki.cgsociety.org/index.php/Ray_Sphere_Intersection
    float a = il_vec3_dot(ray->dir, ray->dir);
    il_vec3 ro = il_vec3_sub(ray->origin, il_positionable_getPosition(&sphere->pos));
    float b = 2 * il_vec3_dot(ray->dir, ro);
    float c = il_vec3_dot(ro, ro) - sphere->radius*sphere->radius;
    float disc = b * b - 4 * a * c;

    if (disc < 0) {
        return false;
    }
    float discsqrt = sqrtf(disc);
    float q = b < 0?
        (-b - discsqrt) / 2.0 :
        (-b + discsqrt) / 2.0;

    float t0 = q / a;
    float t1 = c / q;

    if (t0 > t1) {
        float temp = t0;
        t0 = t1;
        t1 = temp;
    }

    if (t1 < 0) {
        return false;
    }

    out->mat = &sphere->mat;
    float d = out->dist = t0 < 0? t1 : t0;
    float r = sphere->radius;
    out->normal = il_vec3_div(il_vec3_add(il_vec3_mul(ray->dir, il_vec3_new(d,d,d)), ro), il_vec3_new(r,r,r));
    return true;
}

bool test_box(struct box *box, struct ray *ray, struct contact *out)
{
    return false;
}

bool test_light(struct light *light, struct ray *ray, il_vec3 *out)
{
    il_vec3 lp = il_vec3_normal(il_vec3_sub(il_positionable_getPosition(&light->pos), ray->origin));
    float att = il_vec3_dot(ray->dir, lp);
    *out = il_vec3_mul(light->colour, il_vec3_new(att,att,att));
    return att > 0;
}

float clampf(float h, float l, float v)
{
    return v > h? h : v < l? l : v;
}

void set_pixel(il_vec3 *fb, unsigned x, unsigned y, il_vec3 col)
{
    if (x > 800 || y > 600) {
        printf("Coords outside bounds: %u %u\n", x, y);
        abort();
    }
    fb[y*800 + x] = il_vec3_add(fb[y*800 + x], col);
}

// transforms ray into link and spawns 0 or more new rays
void trace(struct ray *ray, struct link *out)
{
    out->id = ray->id;
    out->root = ray->root;
    out->depth = ray->depth;
    out->x = ray->x;
    out->y = ray->y;

    //printf("Tracing ray (%f %f %f) -> (%f %f %f)\n", ray->origin.x, ray->origin.y, ray->origin.z, ray->dir.x, ray->dir.y, ray->dir.z);
    IL_ARRAY(struct contact,) contacts = {0};
    for (unsigned i = 0; i < spheres.length; i++) {
        struct contact c;
        if (test_sphere(&spheres.data[i], ray, &c)) {
            IL_APPEND(contacts, c);
        }
    }
    for (unsigned i = 0; i < boxes.length; i++) {
        struct contact c;
        if (test_box(&boxes.data[i], ray, &c)) {
            IL_APPEND(contacts, c);
        }
    }
    if (contacts.length == 0) {
        il_vec3 col = il_vec3_new(0,0,0);
        for (unsigned i = 0; i < lights.length; i++) {
            il_vec3 out;
            if (test_light(&lights.data[i], ray, &out)) {
                col = il_vec3_add(col, out);
            }
        }
        out->colour = col;
        out->opacity = ray->opacity;
        return;
    }
    float best_dist = INFINITY;
    unsigned best_id = 0;

    for (unsigned i = 0; i < contacts.length; i++) {
        if (contacts.data[i].dist < best_dist) {
            best_dist = contacts.data[i].dist;
            best_id = i;
        }
    }
    struct contact *con = &contacts.data[best_id];
    //printf("Hit! (%f %f %f) -> (%f %f %f), dist: %f, pos: %u, %u\n", ray->origin.x, ray->origin.y, ray->origin.z, ray->dir.x, ray->dir.y, ray->dir.z, contacts.data[best_id].dist, ray->x, ray->y);
    struct ray bounce;
    il_vec3 origin = bounce.origin = il_vec3_add(il_vec3_mul(ray->dir, il_vec3_new(con->dist, con->dist, con->dist)), ray->origin);
    // dir - 2*normal*dot(dir, normal)
    float RdotN = il_vec3_dot(ray->dir, con->normal);
    bounce.dir = il_vec3_sub(
        ray->dir,
        il_vec3_mul(
            il_vec3_new(2,2,2),
        il_vec3_mul(
            con->normal,
            il_vec3_new(RdotN, RdotN, RdotN))));
    bounce.root = ray->root;
    bounce.id = ray_counter++;
    bounce.opacity = il_vec3_mul(ray->opacity, con->mat->specular);
    bounce.depth = ray->depth+1;
    bounce.x = ray->x;
    bounce.y = ray->y;
    if (bounce.depth < max_depth) {
        emit_ray(bounce);
    }
    struct diffuse diffuse;
    diffuse.root = ray->root;
    diffuse.depth = ray->depth+1;
    diffuse.x = ray->x;
    diffuse.y = ray->y;
    diffuse.origin = origin;
    diffuse.opacity = il_vec3_mul(ray->opacity, con->mat->diffuse);
    if (diffuse.depth < max_depth) {
        emit_diffuse(diffuse);
    }

    out->opacity = il_vec3_mul(ray->opacity, con->mat->diffuse);
    out->colour = con->mat->colour;

    IL_FREE(contacts);

    return;
}

const unsigned int num_workers = 4;
struct worker_ctx {
    unsigned num;
    struct ray *rays;
    struct link *links;
    pthread_mutex_t mutex;
    pthread_cond_t input;
} worker_ctxs[num_workers];
pthread_barrier_t worker_barrier;
void *worker(void* ptr)
{
    struct worker_ctx *ctx = ptr;
    pthread_mutex_lock(&ctx->mutex);
    while (1) {
        pthread_cond_wait(&ctx->input, &ctx->mutex);
        for (unsigned i = 0; i < ctx->num; i++) {
            trace(&ctx->rays[i], &ctx->links[i]);
        }
        pthread_barrier_wait(&worker_barrier);
    }
    return NULL;
    pthread_mutex_unlock(&ctx->mutex);
}

il_vec3 *process()
{
    static il_vec3 fb[800*600];
    static il_vec3 out[800*600];
    static struct link links[10000 * num_workers];

    struct ray_chunk *rays_list[4], *cur;
    size_t lens[4] = {0}, len;
    bool have_work = true;
    for (unsigned i = 0; i < num_workers; i++) {
        if (have_work) {
            cur = rays_head;
            rays_head = rays_head->next;
            if (!rays_head) {
                rays_head = rays_tail = ray_chunk_new();
                have_work = false;
            }
            rays_list[i] = cur;
            worker_ctxs[i].num = lens[i] = cur->length;
            worker_ctxs[i].rays = cur->rays;
        } else {
            // we have to signal the other threads regardless because the barrier needs all threads participating
            worker_ctxs[i].num = 0;
        }
        worker_ctxs[i].links = &links[10000 * i];
        pthread_cond_signal(&worker_ctxs[i].input);
    }
    pthread_barrier_wait(&worker_barrier);
    for (unsigned i = 0; i < num_workers; i++) {
        for (unsigned j = 0; j < lens[i]; j++) {
            struct ray *ray = &rays_list[i]->rays[j];
            struct link *link = &links[10000*i + j];
            set_pixel(fb, ray->x, ray->y, il_vec3_mul(link->opacity, link->colour));
        }
        ray_chunk_free(rays_list[i]);
    }

    if (have_work) {
        goto done;
    }
    // if we have no work to do, start hammering on all those diffuse rays
    struct diffuse_chunk *diffuse = diffuse_head;
    diffuse_head = diffuse_head->next;
    if (!diffuse_head) {
        diffuse_head = diffuse_tail = diffuse_chunk_new();
    }
    len = diffuse->length;
    for (unsigned i = 0; i < len; i++) {
        struct diffuse *d = &diffuse->diffuse[i];
        int num_diffuse = 10 - d->depth * 10 / max_depth;
        for (unsigned i = 0; i < num_diffuse; i++) {
            struct ray diffuse;
            // cos-1(2x - 1)
            float a = ((float)rand()/RAND_MAX) * 180;
            float b = acos(((float)rand()/RAND_MAX) * 2 - 1);
            diffuse.dir = il_vec3_new(sin(a), b, cos(a));
            diffuse.origin = d->origin;
            diffuse.id = ray_counter++;
            diffuse.root = d->root;
            diffuse.opacity = il_vec3_mul(d->opacity, il_vec3_new(1.f/num_diffuse, 1.f/num_diffuse, 1.f/num_diffuse));
            diffuse.depth = d->depth + 1;
            diffuse.x = d->x;
            diffuse.y = d->y;
            emit_ray(diffuse);
        }
    }
    diffuse_chunk_free(diffuse);
done:
    ; // it errors if I don't have this here. wtf.
    float sum = 0;
    for (unsigned y = 0; y < 600; y++) {
        for (unsigned x = 0; x < 800; x++) {
            il_vec3 c = fb[y*800 + x];
            sum += il_vec3_len(c);
        }
    }
    float avg = 800*600 / sum;
    for (unsigned y = 0; y < 600; y++) {
        for (unsigned x = 0; x < 800; x++) {
            out[y*800+x] = il_vec3_mul(fb[y*800+x], il_vec3_new(avg,avg,avg));
        }
    }

    return out;
}

void setupScene()
{
    rays_head = rays_tail = ray_chunk_new();
    diffuse_head = diffuse_tail = diffuse_chunk_new();
    world = il_world_new(1);

    unsigned int seed = 0xDEADBEEF;
    for (unsigned i = 0; i < 10; i++) {
        struct sphere test;
        test.pos = il_positionable_new(world);
        il_positionable_setPosition(&test.pos, il_vec3_new(
            ((float)rand_r(&seed) / RAND_MAX) * 4 - 2,
            ((float)rand_r(&seed) / RAND_MAX) * 4 - 2,
            ((float)rand_r(&seed) / RAND_MAX) * 4 + 7));
        test.radius = 1;
        test.mat.colour = il_vec3_new(
            (float)rand_r(&seed) / RAND_MAX,
            (float)rand_r(&seed) / RAND_MAX,
            (float)rand_r(&seed) / RAND_MAX
        );
        test.mat.diffuse = il_vec3_mul(test.mat.colour, il_vec3_new(.5, .5, .5));
        test.mat.specular = test.mat.colour;
        test.mat.specular_co = 96;
        IL_APPEND(spheres, test);
    }

    struct light l;
    l.pos = il_positionable_new(world);
    il_positionable_setPosition(&l.pos, il_vec3_new(-3, 3, 7));
    l.radius = 10;
    l.colour = il_vec3_new(1.0, 1.0, 0.0);
    IL_APPEND(lights, l);

    int w = 800*4, h = 600*4;
    for (int y = 0; y < h; y++) {
        for (int x = 0; x < w; x++) {
            struct link link;
            link.id = ray_counter++;
            link.opacity = il_vec3_new(1,1,1);
            link.root = link.id; // rooted to itself
            //root_links[y * 800 + x] = link;
            struct ray ray;
            ray.id = ray_counter++;
            ray.root = link.root;
            ray.origin = il_vec3_new(0,0,0);
            ray.opacity = il_vec3_new(.125,.125,.125);
            ray.depth = 0;
            ray.x = x/4;
            ray.y = y/4;
            // projection
            ray.dir = il_vec3_normal(il_vec3_new((x - w/2) / (float)w, (y - h/2) / (float)h, 1));
            emit_ray(ray);
        }
    }
}

const char vert_source[] =
"#version 140\n"
"in vec2 in_Position;\n"
"out vec2 texcoord;\n"
"void main()\n"
"{\n"
"   gl_Position = vec4(in_Position, 0, 0);\n"
"   texcoord = in_Position;\n"
"}\n";

const char frag_source[] =
"#version 140\n"
"in vec2 texcoord;\n"
"uniform sampler2D tex;\n"
"out vec4 out_Color;\n"
"void main()\n"
"{\n"
"   out_Color = texture(tex, texcoord);\n"
"}\n";

int main(int argc, char **argv)
{
    (void)argc,(void)argv;
    //il_opts opts = il_opts_parse(argc, argv);
    il_add_module_path("../IntenseLogic/build");
    il_load_ilutil();
    il_load_ilmath();
    il_load_ilcommon();
    il_load_ilgraphics();
    pthread_t workers[num_workers];
    pthread_mutex_init(&diffuse_mutex, NULL);
    pthread_mutex_init(&ray_mutex, NULL);
    pthread_barrier_init(&worker_barrier, NULL, num_workers + 1);
    for (unsigned i = 0; i < num_workers; i++) {
        pthread_create(&workers[i], NULL, worker, &worker_ctxs[i]);
        pthread_cond_init(&worker_ctxs[i].input, NULL);
        pthread_mutex_init(&worker_ctxs[i].mutex, NULL);
    }
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_MAJOR_VERSION, 2);
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_MINOR_VERSION, 1);
    SDL_GL_SetAttribute(SDL_GL_RED_SIZE, 8);
    SDL_GL_SetAttribute(SDL_GL_GREEN_SIZE, 8);
    SDL_GL_SetAttribute(SDL_GL_BLUE_SIZE, 8);
    /*SDL_GL_SetAttribute(SDL_GL_CONTEXT_FLAGS, SDL_GL_CONTEXT_DEBUG_FLAG | SDL_GL_CONTEXT_FORWARD_COMPATIBLE_FLAG);
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_PROFILE_MASK, SDL_GL_CONTEXT_PROFILE_CORE);*/
    SDL_Window *window;
    if (!(window = SDL_CreateWindow("Ray Tracer", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED, 800, 600, SDL_WINDOW_OPENGL | SDL_WINDOW_RESIZABLE))) {
        printf("Couldn't open window\n");
        return 1;
    }
    SDL_GLContext context;
    if (!(context = SDL_GL_CreateContext(window))) {
        printf("Couldn't create context\n");
        return 1;
    }
    glewExperimental = GL_TRUE;
    if (glewInit() != GLEW_OK) {
        printf("Couldn't init GLEW\n");
        return 1;
    }

    ilG_testError("glew");
    GLuint tex, vbo, vao, prog, vert, frag;
    glGenTextures(1, &tex);
    /*glGenBuffers(1, &vbo);
    glGenVertexArrays(1, &vao);*/

    glEnable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D, tex);
    glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    /*glBindBuffer(GL_ARRAY_BUFFER, vbo);
    float data[] = {
        -1,-1,  -1,+1,  +1,-1,
        +1,-1,  -1,+1,  +1,+1
    };
    glBufferData(GL_ARRAY_BUFFER, sizeof(data), data, GL_STATIC_DRAW);
    glBindVertexArray(vao);
    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, NULL);
    glEnableVertexAttribArray(0);
    ilG_testError("Setting up vbo");
    prog = glCreateProgram();
    vert = ilG_makeShader(GL_VERTEX_SHADER, il_l(vert_source));
    frag = ilG_makeShader(GL_FRAGMENT_SHADER, il_l(frag_source));
    glAttachShader(prog, vert);
    glAttachShader(prog, frag);
    glBindAttribLocation(prog, 0, "in_Position");
    glBindFragDataLocation(prog, 0, "out_Color");
    ilG_testError("Setting up shader");
    if (ilG_linkProgram(prog)) {
        printf("Couldn't link shader\n");
        return 1;
    }
    glUseProgram(prog);
    glUniform1i(glGetUniformLocation(prog, "tex"), GL_TEXTURE0);
    ilG_testError("glUniform");*/

    setupScene();

    struct timeval start, end;
    unsigned int frame = 0;
    gettimeofday(&start, NULL);
    //SDL_GL_SetSwapInterval(0);
    int processing = 1;
    unsigned max_memory = 0;
    while (1) {
        frame++;
        SDL_Event ev;
        while (SDL_PollEvent(&ev)) {
            switch (ev.type) {
            case SDL_QUIT: 
                {
                    gettimeofday(&end, NULL);
                    float sec = (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1000000.f;
                    printf("%u frames over %f seconds: %f fps\n", frame, sec, frame / sec);
                    printf("%.4fmb peak memory usage\n", max_memory / 1000000.f);
                    for (unsigned i = 0; i < num_workers; i++) {
                        pthread_cancel(workers[i]);
                    }
                    return 0;
                }
            case SDL_KEYDOWN:
                processing = !processing;
                break;
            default:
                break;
            }
        }
        glClearColor(0,0,0,0);
        //glClear(GL_COLOR_BUFFER_BIT);
        //glDrawArrays(GL_TRIANGLES, 0, 6);
        if (processing) {
            unsigned bytes = 0, count = 0;
            struct ray_chunk *rcur = rays_head;
            while (rcur) {
                bytes += sizeof(struct ray_chunk);
                count += rcur->length;
                rcur = rcur->next;
            }
            rcur = rays_free;
            while (rcur) {
                bytes += sizeof(struct ray_chunk);
                rcur = rcur->next;
            }
            printf("%09u rays, ", count);
            struct diffuse_chunk *dcur = diffuse_head;
            count = 0;
            while (dcur) {
                bytes += sizeof(struct diffuse_chunk);
                count += dcur->length;
                dcur = dcur->next;
            }
            dcur = diffuse_free;
            while (dcur) {
                bytes += sizeof(struct diffuse_chunk);
                dcur = dcur->next;
            }
            printf("%09u diffuse, ", count);
            count = 800*600 + links.length;
            bytes += count * sizeof(struct link);
            if (bytes > max_memory) {
                max_memory = bytes;
            }
            printf("%09u links, %.4fmb of memory\n", count, bytes / 1000000.f);

            il_vec3 *fb = process();
            glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, 800, 600, 0, GL_RGBA, GL_FLOAT, fb);
            glBegin(GL_TRIANGLE_FAN);
            glVertex2f  (-1, -1);
            glTexCoord2f( 0,  0);
            glVertex2f  (-1, +1);
            glTexCoord2f( 0,  1);
            glVertex2f  (+1, +1);
            glTexCoord2f( 1,  1);
            glVertex2f  (+1, -1);
            glTexCoord2f( 1,  0);
            glEnd();
        }

        ilG_testError("glDrawArrays");
        SDL_GL_SwapWindow(window);
    }
    return 0;
}