6
\$\begingroup\$

In the past I had followed the Ray Tracing in a Weekend books using C++ because that was what the book used. However, recently I started thinking about how hard it would be to implement using C and whether I could do it. As a result I set a goal for implementing it using only C and I was able to mostly accomplish that, except for 1 file written in C++.

Compile with GCC, with extensions enabled (e.g. gcc -std=gnu11).

ray_trace.c

/* standard headers */
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <float.h>
#include <time.h>

/* SDL2 headers */
#include <SDL2/SDL.h>

/* common project headers */
#include "types.h"
#include "vec.h"
#include "ray_trace.h"

/* constants */
#define ASPECT_RATIO (1.0f)
enum { WIDTH = 600, HEIGHT = (i32)(WIDTH / ASPECT_RATIO) };
enum { SAMPLES_PER_PIXEL = 1000 };

/* main project headers */
#include "ray.h"
#include "camera.h"
#include "texture.h"
#include "aabb.h"
#include "hittable.h"
#include "material.h"
#include "hittable_list.h"
#include "sphere.h"
#include "moving_sphere.h"
#include "aarect.h"
#include "box.h"
#include "constanst_medium.h"
#include "bvh.h"

internal float3 ray_color(ray const *r, float3 background,  hittable const *world, i32 depth)
{
    hit_record rec;

    /* if we have exceeded the ray bounce limit, no more light is gathered. */
    if(depth <= 0)
        return 0.0f;
    
    /* if we hit something */
    if (!world->vtable->hit(world, r, 0.001f, FLT_MAX, &rec)) 
        return background;
    
    ray scattered;
    float3 emitted = rec.mat_ptr->vtable->emitted(rec.mat_ptr, rec.uv, rec.p);
    float3 attenuation;

    if(!rec.mat_ptr->vtable->scatter(rec.mat_ptr, r, &rec, &attenuation, &scattered))
        return emitted;

    return emitted + attenuation * ray_color(&scattered, background, world, depth - 1);
}

internal hittable *random_scene(void)
{
    hittable *world = new(hittable_list);

    /* preallocate memory for the world */
    hittable_list_reserve(world, 22 * 22);

    texture *checker = new_checker_texture((float3){0.2f, 0.3f, 0.1f}, 0.9f);
    material *material_ground = new_lambertian(checker);
    hittable_list_add(world, new(sphere, (float3){0, -1000, -1}, 1000, material_ground));

    for(i32 a = -11; a < 11; ++a)
    {
        for(i32 b = -11; b < 11; ++b)
        {
            float choose_mat = randomf();
            float3 center = {a + 0.9f * randomf(), 0.2f, b + 0.9f * randomf()};

            if(length(center - (float3){4, 0.2f, 0}) > 0.9f)
            {
                /* diffuse */
                if(choose_mat < 0.8f)
                {
                    float3 albedo = randomf3() * randomf3();
                    material *sphere_material = new_lambertian(albedo);
                    hittable_list_add(world, new(sphere, center, 0.2f, sphere_material));
                }
                /* metal */
                else if(choose_mat < 0.95f)
                {
                    float3 albedo = randomf3(0.5f, 1.0f);
                    float fuzz = randomf(0.0f, 0.5f);
                    material *sphere_material = new_metal(albedo, fuzz);
                    hittable_list_add(world, new(sphere, center, 0.2f, sphere_material));
                }
                /* glass */
                else
                {
                    material *sphere_material = new_dielectric(1.5f);
                    hittable_list_add(world, new(sphere, center, 0.2f, sphere_material));
                }
            }
        }
    }
    material *material1 = new_dielectric(1.5f);
    hittable_list_add(world, new(sphere, (float3){0, 1, 0}, 1.0f, material1));

    material *material2 = new_lambertian((float3){0.4f, 0.2f, 0.1f});
    hittable_list_add(world, new(sphere, (float3){-4, 1, 0}, 1.0f, material2));

    material *material3 = new_metal((float3){0.7f, 0.6f, 0.6f}, 0.0f);
    hittable_list_add(world, new(sphere, (float3){4, 1, 0}, 1.0f, material3));
    world = new_bvh_node(world, 0, 0);
        
    return world;
}

internal hittable *two_spheres(void)
{
    hittable *objects = new(hittable_list);
    texture *checker = new_checker_texture((float3){0.2f, 0.3f, 0.1f}, 0.9f);

    hittable_list_add(objects, new(sphere, (float3){0, -10, 0}, 10, new_lambertian(checker)));
    hittable_list_add(objects, new(sphere, (float3){0, 10, 0}, 10, new_lambertian(checker)));

    return objects;
}

internal hittable *two_perlin_spheres(void)
{
    hittable *objects = new(hittable_list);

    texture *pertext = new_noise_texture(4);
    hittable_list_add(objects, new(sphere, (float3){0, -1000, 0}, 1000, new_lambertian(pertext)));
    hittable_list_add(objects, new(sphere, (float3){0, 2, 0}, 2, new_lambertian(pertext)));

    return objects;
}

internal hittable *simple_light(void)
{
    hittable *objects = new(hittable_list);

    texture *pertext = new_noise_texture(4);
    hittable_list_add(objects, new(sphere, (float3){0, -1000, 0}, 1000, new_lambertian(pertext)));
    hittable_list_add(objects, new(sphere, (float3){0, 2, 0}, 2, new_dielectric(1.4f)));
    material *difflight = new_diffuse_light(4);
    hittable_list_add(objects, new_xy_rect(3, 5, 1, 3, -2, difflight));
    hittable_list_add(objects, new(sphere, (float3){0, 5, 0}, 1.5f, difflight));


    return objects;
}

internal hittable *cornell_box(void)
{
    hittable *objects = new(hittable_list);

    material *red = new_lambertian((float3){0.65f, 0.05f, 0.05f});
    material *white = new_lambertian(.73f);
    material *green = new_lambertian((float3){0.12f, 0.45f, 0.15f});
    material *light = new_diffuse_light((float3){15, 15, 15});

    /* add walls and light */
    hittable_list_add(objects, new_yz_rect(0, 555, 0, 555, 555, green));
    hittable_list_add(objects, new_yz_rect(0, 555, 0, 555,  0, red));
    hittable_list_add(objects, new_xz_rect(113, 443, 127, 432, 554, light));
    hittable_list_add(objects, new_xz_rect(0, 555, 0, 555,  0, white));
    hittable_list_add(objects, new_xz_rect(0, 555, 0, 555, 555, white));
    hittable_list_add(objects, new_xy_rect(0, 555, 0, 555, 555, white));

    /* add box1 */
    hittable *box1 = new(box, 0, (float3){165, 330, 165}, white);
    box1 = new(rotate_y, box1, 15);
    box1 = new(translate, box1, (float3){265, 0, 295});
    hittable_list_add(objects, box1);
    
    /* add box2 */
    hittable *box2 = new(box, 0, 165, white);
    box2 = new(rotate_y, box2, -18);
    box2 = new(translate, box2, (float3){130, 0, 65});
    hittable_list_add(objects, box2);
    
    
    return objects;
}

internal hittable *cornell_smoke(void)
{
    hittable *objects = new(hittable_list);

    material *red = new_lambertian((float3){0.65f, 0.05f, 0.05f});
    material *white = new_lambertian(.73f);
    material *green = new_lambertian((float3){0.12f, 0.45f, 0.15f});
    material *light = new_diffuse_light((float3){7, 7, 7});

    /* add walls and light */
    hittable_list_add(objects, new_yz_rect(0, 555, 0, 555, 555, green));
    hittable_list_add(objects, new_yz_rect(0, 555, 0, 555,  0, red));
    hittable_list_add(objects, new_xz_rect(113, 443, 127, 432, 554, light));
    hittable_list_add(objects, new_xz_rect(0, 555, 0, 555,  0, white));
    hittable_list_add(objects, new_xz_rect(0, 555, 0, 555, 555, white));
    hittable_list_add(objects, new_xy_rect(0, 555, 0, 555, 555, white));

    /* create box1 */
    hittable *box1 = new(box, 0, (float3){165, 330, 165}, white);
    box1 = new(rotate_y, box1, 15);
    box1 = new(translate, box1, (float3){265, 0, 295});

    /* create box2 */
    hittable *box2 = new(box, 0, 165, white);
    box2 = new(rotate_y, box2, -18);
    box2 = new(translate, box2, (float3){130, 0, 65});

    /* add boxes */
    hittable_list_add(objects, new(constant_medium, box1, 0.01f, 0.0f));
    hittable_list_add(objects, new(constant_medium, box2, 0.01f, 1.0f));

    return objects;
}

internal hittable *final_scene(void)
{
    hittable *boxes1 = new(hittable_list);
    material *ground = new_lambertian((float3){0.48f, 0.83f, 0.53f});

    i32 const boxes_per_side = 20;
    hittable_list_reserve(boxes1, boxes_per_side * boxes_per_side);
    for(i32 i = 0; i < boxes_per_side; ++i)
    {
        for (i32 j = 0; j < boxes_per_side; ++j)
        {
            float w = 100.0f;
            float x0 = -1000.0f + i * w;
            float z0 = -1000.0f + j * w;
            float y0 = 0.0f;
            float x1 = x0 + w;
            float y1 = randomf(1, 101);
            float z1 = z0 + w;

            hittable_list_add(boxes1, new(box, (float3){x0, y0, z0}, (float3){x1, y1, z1}, ground));
        }
    }

    hittable *objects = new(hittable_list, new_bvh_node(boxes1, 0, 1));

    material *light = new_diffuse_light((float3){7, 7, 7});
    hittable_list_add(objects, new_xz_rect(123, 423, 147, 412, 554, light));

    float3 center1 = {400, 400, 200};
    float3 center2 = center1 + (float3){30, 0, 0};
    material *moving_sphere_material = new_lambertian((float3){0.7f, 0.3f, 0.1f});
    hittable_list_add(objects, new_moving_sphere(center1, center2, 0, 1, 50, moving_sphere_material));

    hittable_list_add(objects, new(sphere, (float3){260, 150, 45}, 50, new_dielectric(1.5f)));
    hittable_list_add(objects, new(sphere, (float3){0, 150, 145}, 50, new_metal((float3){0.8f, 0.8f, 0.9f}, 10.0f)));

    hittable *boundary = new(sphere, (float3){360, 150, 145}, 70, new_dielectric(1.5f));
    hittable_list_add(objects, boundary);
    hittable_list_add(objects, new(constant_medium, boundary, 0.2f, (float3){0.2f, 0.4f, 0.9f}));
    boundary = new(sphere, 0, 5000, new_dielectric(1.5f));
    hittable_list_add(objects, new(constant_medium, boundary, 0.0001f, (float3){1, 1, 1}));

    material *emat = new_lambertian(new_solid_color((float3){0, 0, 0.75f}));
    hittable_list_add(objects, new(sphere, (float3){400, 200, 400}, 100, emat));
    texture *pertext = new_noise_texture(0.1f);
    hittable_list_add(objects, new(sphere, (float3){220, 280, 300}, 80, new_lambertian(pertext)));

    hittable *boxes2 = new(hittable_list);
    material *white = new_lambertian(.73f);
    i32 const ns = 1000;
    hittable_list_reserve(boxes2, ns);
    for(i32 j = 0; j < ns; ++j)
    {
        hittable_list_add(boxes2, new(sphere, randomf3(0, 165), 10, white));
    }

    hittable_list_add(objects, new(translate, new(rotate_y, new_bvh_node(boxes2, 0.0f, 1.0f), 15), (float3){-100, 270, 395}));

    return objects;
}

internal void generate_image(uint32_t *image)
{
    hittable *world;
    float3 lookfrom;
    float3 lookat;
    float aperture;
    float vfov;
    float3 background;
    
    /* choose the scene */
    switch (7)
    {
        case 1:
            world = random_scene();
            background = (float3){0.70f, 0.80f, 1};
            lookfrom = (float3){13, 2, 3};
            lookat = (float3){0, 0, 0};
            vfov = 20.0f;
            aperture = 0.0f;
            break;
        case 2:
            world = two_spheres();
            background = (float3){0.70f, 0.80f, 1};
            lookfrom = (float3){13, 2, 3};
            lookat = 0;
            vfov = 20.0f;
            aperture = 0.0f;
            break;
        case 3:
            world = two_perlin_spheres();
            background = (float3){0.70f, 0.80f, 1};
            lookfrom = (float3){13, 2, 3};
            lookat = 0;
            vfov = 20.0f;
            aperture = 0.0f;
            break;
        case 4:
            world = simple_light();
            background = 0;
            lookfrom = (float3){26, 3, 6};
            lookat = (float3){0, 2, 0};
            vfov = 20.0f;
            aperture = 0.0f;
            break;
        case 5:
            world = cornell_box();
            background = 0;
            lookfrom = (float3){278, 278, -800};
            lookat = (float3){278, 278, 0};
            vfov = 40.0f;
            aperture = 0.0f;
            break;
        case 6:
            world = cornell_smoke();
            background = 0;
            lookfrom = (float3){278, 278, -800};
            lookat = (float3){278, 278, 0};
            vfov = 40.0f;
            aperture = 0.0f;
            break;
        default:
        case 7:
            world = final_scene();
            background = 0;
            lookfrom = (float3){478, 278, -600};
            lookat = (float3){278, 278, 0};
            vfov = 40.0f;
            aperture = 0.0f;
            break;
    }

    /* setup camera */
    float3 vup = (float3){0, 1, 0};
    float dist_to_focus = 10;
    camera cam = make_camera(lookfrom, lookat, vup, vfov, ASPECT_RATIO, aperture, dist_to_focus, 0, 1.0f);
        
    
    /* render */
    #pragma omp parallel for
    for(i32 pix = 0; pix < HEIGHT * WIDTH; ++pix) /* pixel index */
    {   
        float3 pixel_color = {0};
        for(i32 s = 0; s < SAMPLES_PER_PIXEL; ++s)
        {
            /* get pixel coords */
            int32_t i = pix % WIDTH;
            int32_t j = (HEIGHT - 1) - (pix / WIDTH);
        
            float u = (i + randomf())/(float)(WIDTH - 1);
            float v = (j + randomf())/(float)(HEIGHT - 1);
            ray r = camera_get_ray(&cam, u, v);

            /* add to the samples */
            {
                float3 result_color = ray_color(&r, background, world, 50);
                pixel_color += result_color;    
            }
        }

        write_pixel(&image[pix], pixel_color, SAMPLES_PER_PIXEL);
    }
}

int main(int argc, char *argv[])
{
    /* we don't need these
     * but SDL2 needs main to always have the signature of int main(int, char*[]);
     */
    (void)argc;
    (void)argv;

    /* setup SDL2 */
    if(SDL_Init(SDL_INIT_VIDEO) < 0)
    {
        fprintf(stderr, "Error: with SDL2 initialization \"%s\"\n", SDL_GetError());
        return -1;
    }

    SDL_Window *window = SDL_CreateWindow("", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED, WIDTH, HEIGHT, SDL_WINDOW_OPENGL);
    if(NULL == window)
    {
        fprintf(stderr, "Error: with SDL2 window creation \"%s\"\n", SDL_GetError());
        return -2;
    }

    SDL_Renderer *renderer = SDL_CreateRenderer(window, -1, SDL_RENDERER_ACCELERATED);
    if(NULL == renderer)
    {
        fprintf(stderr, "Error: with SDL2 renderer creation \"%s\"\n", SDL_GetError());
        return -3;  
    }

    SDL_Texture *texture = SDL_CreateTexture(renderer, SDL_PIXELFORMAT_RGBA8888, SDL_TEXTUREACCESS_STATIC, WIDTH, HEIGHT);
    if(NULL == texture)
    {
        fprintf(stderr, "Error: with SDL2 texture creation \"%s\"\n", SDL_GetError());
        return -4;  
    }

    /* allocate memory for an image */
    uint32_t *image = malloc(sizeof(uint32_t) * WIDTH * HEIGHT);

    /* generate the image and time how long it took */
    clock_t start_time = clock();
    generate_image(image);      
    clock_t end_time = clock();

    printf("it took %f seconds", (float)(end_time - start_time) / (float)CLOCKS_PER_SEC);

    /*  main loop */
    bool quit = false;
    while(!quit)
    {
        /* check for a quit event */
        {
            SDL_Event event = { 0 };
            if(SDL_PollEvent(&event))
                switch (event.type)
                {
                    case SDL_QUIT:
                    {
                        quit = true;
                    } break;
                }
        }
        
        /* draw texture */
        {
            SDL_UpdateTexture(texture, NULL, image, sizeof(u32) * WIDTH);
            SDL_RenderCopy(renderer, texture, NULL, NULL);
            SDL_RenderPresent(renderer);
        }
    }

    /* success */
    return 0;
}

ray_trace.h

#pragma once
#define PI (3.141592653589f)
#define new(type, ...) ({type *$$$$ = malloc(sizeof(type)); *$$$$ = make_##type(__VA_ARGS__); (void*)$$$$;})

float random_float(void);

internal inline void write_pixel(u32 *pixel, float3 color, i32 sample_count)
{
    /* average the samples and do gamma correction */
    {
        color /= sample_count;
        color[0] = sqrt(color[0]);
        color[1] = sqrt(color[1]);
        color[2] = sqrt(color[2]);
    }
    *pixel = ((u32)(fminf(color[0] * 255.0f + 0.5555f, 255.999f)) << 24) | ((u32)(fminf(color[1] * 255.0f + 0.5555f, 255.999f)) << 16) | ((u32)(fminf(color[2] * 255.0f + 0.5555f, 255.999f)) << 8);
}

internal inline size_t index_at(size_t x, size_t y, size_t width)
{
    return y  * width + x; 
}

internal inline float degrees_to_radians(float degrees) 
{
    return degrees * PI / 180.0f;
}

internal inline float clamp(float x, float min, float max) {
    if (x < min) return min;
    if (x > max) return max;
    return x;
}

/* random number genreration */
__attribute__((overloadable)) internal inline float randomf(void)
{
    return random_float();
}

__attribute__((overloadable)) internal inline float randomf(float min, float max)
{
    return min + (max - min) * randomf();
}

__attribute__((overloadable)) internal inline float2 randomf2(void)
{
    return (float2){randomf(), randomf()};
}

__attribute__((overloadable)) internal inline float3 randomf3(void)
{
    return (float3){randomf(), randomf(), randomf()};
}

__attribute__((overloadable)) internal inline float4 randomf4(void)
{
    return (float4){randomf(), randomf(), randomf(), randomf()};
}

__attribute__((overloadable)) internal inline float2 randomf2(float min, float max)
{
    return (float2){randomf(min, max), randomf(min, max)};
}

__attribute__((overloadable)) internal inline float3 randomf3(float min, float max)
{
    return (float3){randomf(min, max), randomf(min, max), randomf(min, max)};
}

__attribute__((overloadable)) internal inline float4 randomf4(float min, float max)
{
    return (float4){randomf(min, max), randomf(min, max), randomf(min, max), randomf(min, max)};
}

internal inline float3 random_in_unit_sphere(void)
{
    for(;;)
    {
        float3 p = randomf3(-1, 1);
        if(length_sqr(p) >= 1) continue;
        return p;
    }
}

internal inline float3 random_unit_vector(void) {
    float a = randomf(0, 2 * PI);
    float z = randomf(-1, 1);
    float r = sqrtf(1 - z*z);
    return (float3){r * cosf(a), r * sinf(a), z};
}

internal inline float3 random_in_unit_disk(void)
{
    for(;;)
    {
        float3 p = {randomf(-1, 1), randomf(-1, 1), 0};
        if(length_sqr(p) >= 1) continue;
        return p;
    }
}

vec.h

#pragma once
#include <math.h>

/* types */
typedef float float2 __attribute__((ext_vector_type(2)));
typedef float float3 __attribute__((ext_vector_type(3)));
typedef float float4 __attribute__((ext_vector_type(4)));

/* length squared */
__attribute__((overloadable)) internal inline float length_sqr(float2 v)
{
    return v.x * v.x + v.y * v.y; 
}

__attribute__((overloadable)) internal inline float length_sqr(float3 v)
{
    return v.x * v.x + v.y * v.y + v.z * v.z; 
}

__attribute__((overloadable)) internal inline float length_sqr(float4 v)
{
    return v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w; 
}

/* length */
__attribute__((overloadable)) internal inline float length(float2 v)
{
    return sqrtf(length_sqr(v)); 
}

__attribute__((overloadable)) internal inline float length(float3 v)
{
    return sqrtf(length_sqr(v)); 
}

__attribute__((overloadable)) internal inline float length(float4 v)
{
    return sqrtf(length_sqr(v)); 
}

/* dot product */
__attribute__((overloadable)) internal inline float dot(float2 u, float2 v)
{
    return u[0] * v[0] + u[1] * v[1]; 
}

__attribute__((overloadable)) internal inline float dot(float3 u, float3 v)
{
    return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; 
}

__attribute__((overloadable)) internal inline float dot(float4 u, float4 v)
{
    return u[0] * v[0] + u[1] * v[1] + u[2] * v[2] + v[3] * v[3]; 
}

/* cross product */
__attribute__((overloadable)) internal inline float3 cross(float3 u, float3 v)
{
    return (float3){u[1] * v[2] - u[2] * v[1],
                    u[2] * v[0] - u[0] * v[2],
                    u[0] * v[1] - u[1] * v[0]};
}

/* normalization */
__attribute__((overloadable)) internal inline float2 normalize(float2 v)
{
    return v / length(v);
}

__attribute__((overloadable)) internal inline float3 normalize(float3 v)
{
    return v / length(v);
}

__attribute__((overloadable)) internal inline float4 normalize(float4 v)
{
    return v / length(v);
}

/* reflection */
internal inline float3 reflect(float3 v, float3 n)
{
    return v - 2 * dot(v, n) * n;
}

/* refraction */
internal inline float3 refract(float3 uv, float3 n, float etai_over_etat)
{
    float cos_theta = dot(-uv, n);
    float3 r_out_perp = etai_over_etat * (uv + cos_theta * n);
    float3 r_out_parallel = -sqrtf(fabsf(1.0f - length_sqr(r_out_perp))) * n;
    return r_out_parallel + r_out_perp;
}

types.h

#pragma once
#include <stdint.h>

/* integral types */
typedef int8_t i8;
typedef int16_t i16;
typedef int32_t i32;
typedef int64_t i64;

/* unsigned integral types */
typedef uint8_t u8;
typedef uint16_t u16;
typedef uint32_t u32;
typedef uint64_t u64;

/*  floating point types */
typedef float f32;
typedef double f64;

/* static has a lot of uses */
#define internal static
#define auto __auto_type

ray.h

#pragma once
typedef struct ray
{
    float3 pos;
    float3 dir;
    float time;
} ray;

internal inline float3 ray_at(ray const *this, float t)
{
    return this->pos + t * this->dir;
}

camera.h

#pragma once

typedef struct camera
{
    float3 origin;
    float3 horizontal;
    float3 vertical;
    float3 lower_left_corner;
    float3 u, v, w;
    float lens_radius;
    float time0, time1;
} camera;

internal inline ray camera_get_ray(camera const *this, float s, float t)
{
    float3 rd = this->lens_radius * random_in_unit_disk();
    float3 offset = this->u * rd.x + this->v * rd.y;

    return (ray){
        this->origin + offset,
        this->lower_left_corner + s * this->horizontal + t * this->vertical - this->origin - offset,
        randomf(this->time0, this->time1)
    };
}

internal camera make_camera(float3 lookfrom, float3 lookat, float3 vup, float fov, float aspect_ratio, float aperture, float focus_dist, float t0, float t1)
{
    float theta = degrees_to_radians(fov);
    float h = tanf(theta / 2);
    float viewport_height = 2.0f * h;
    float viewport_width = aspect_ratio * viewport_height;

    float3 w = normalize(lookfrom - lookat);
    float3 u = normalize(cross(vup, w));
    float3 v = cross(w, u);
    camera result = {
        .w = w,
        .u = u,
        .v = v,
        .origin = lookfrom,
        .horizontal = focus_dist * viewport_width * u,
        .vertical= focus_dist * viewport_height * v,
        .lower_left_corner = result.origin - result.horizontal / 2 - result.vertical / 2 - w * focus_dist,
        .lens_radius = aperture / 2,
        .time0 = t0,
        .time1 = t1
    };

    return result;      
}

perlin.h

#pragma once

enum { PERLIN_POINT_COUNT = 256 };

typedef struct perlin
{
    float3 *ranvec;
    i32 *perm_x;
    i32 *perm_y;
    i32 *perm_z;
} perlin;

internal void perlin_permute(i32 *p, i32 n)
{
    for(i32 i = n - 1; i > 0; --i)
    {
        i32 target = (i32)randomf(0, i + 1);
        i32 temp = p[i];
        p[i] = p[target];
        p[target] = temp;
    }   
}

internal i32 *perlin_generate_perm(void)
{
    i32 *p = malloc(sizeof(i32) * PERLIN_POINT_COUNT);
    for(i32 i = 0; i < PERLIN_POINT_COUNT; ++i)
        p[i] = i;
        
    perlin_permute(p, PERLIN_POINT_COUNT);

    return p;
}

internal float perlin_interp(float3 c[2][2][2], float u, float v, float w)
{
    float uu = u * u * (3 - 2 * u);
    float vv = v * v * (3 - 2 * v);
    float ww = w * w * (3 - 2 * w);
    float result = 0.0f;

    for(i32 i = 0; i < 2; ++i)
        for(i32 j = 0; j < 2; ++j)
            for(i32 k = 0; k < 2; ++k)
            {
                float3 weight_v = { u - i, v - j, w - k};
                result += ((i * uu) + (1 - i) * (1 - uu))
                        * ((j * vv) + (1 - j) * (1 - vv))
                        * ((k * ww) + (1 - k) * (1 - ww))
                        * dot(c[i][j][k], weight_v);
            }
    return result;
}

internal float perlin_noise(perlin const *this, float3 p)
{
    float u = p.x - floorf(p.x);
    float v = p.y - floorf(p.y);
    float w = p.z - floorf(p.z);

    i32 i = floorf(p.x);
    i32 j = floorf(p.y);
    i32 k = floorf(p.z);
    float3 c[2][2][2];

    for(i32 di = 0; di < 2; ++di)
        for(i32 dj = 0; dj < 2; ++dj)
            for(i32 dk = 0; dk < 2; ++dk)
            {
                c[di][dj][dk] =
                    this->ranvec[
                        this->perm_x[(i + di) & 255] ^
                        this->perm_y[(j + dj) & 255] ^
                        this->perm_z[(k + dk) & 255]
                    ];
            }

    return perlin_interp(c, u, v, w);
}

internal float perlin_turb(perlin const *this, float3 p, i32 depth)
{
    float result = 0;
    float3 temp_p = p;
    float weight = 1;

    for(i32 i = 0; i < depth; ++i)
    {
        result += weight * perlin_noise(this, temp_p);
        weight *= 0.5f;
        temp_p *= 2;
    }

    return  fabsf(result);
}

internal perlin make_perlin(void)
{
    perlin result = { .ranvec = malloc(sizeof(float3) * PERLIN_POINT_COUNT) };
    for(i32 i = 0; i < PERLIN_POINT_COUNT; ++i)
    {
        result.ranvec[i] = normalize(randomf3(-1, 1));
    }
    
    result.perm_x = perlin_generate_perm();
    result.perm_y = perlin_generate_perm();
    result.perm_z = perlin_generate_perm();

    return result;
}

texture.h

#pragma once
#include "perlin.h"

struct texture;

struct texture_vtable
{
    float3(*value)(struct texture const *this, float2 uv, float3 p);
};

typedef struct texture
{
    struct texture_vtable *vtable;  
} texture;

typedef struct solid_color
{
    texture parent;
    float3 color_value;
} solid_color;

internal float3 solid_color_value(solid_color const *this, float2 uv, float3 p)
{
    return this->color_value;
}

internal texture *new_solid_color(float3 c)
{
    solid_color *result = malloc(sizeof(solid_color));
    result->color_value = c;

    /* setup vtable */
    static struct texture_vtable solid_color_vtable = { .value = (void*)solid_color_value };
    result->parent.vtable = &solid_color_vtable;

    return (void*)result;
}

typedef struct checker_texture
{
    texture parent;
    texture *odd;
    texture *even;
} checker_texture;

internal float3 checker_texture_value(checker_texture const *this, float2 uv, float3 p)
{
    float sines = sinf(10 * p.x) * sinf(10 * p.y) * sinf(10 * p.z);
    if(sines < 0)
        return this->odd->vtable->value(this->odd, uv, p);
    else
        return this->even->vtable->value(this->even, uv, p);
}

__attribute__((overloadable)) internal texture *new_checker_texture(texture *t0, texture *t1)
{
    checker_texture *result = malloc(sizeof(checker_texture));
    result->even = t0;
    result->odd = t1;

    /* setup vtable */
    static struct texture_vtable checker_texture_vtable = { .value  = (void*)&checker_texture_value };
    result->parent.vtable = &checker_texture_vtable;

    return (void*)result;
}

__attribute__((overloadable)) internal texture *new_checker_texture(float3 c1, float3 c2)
{
    return new_checker_texture(new_solid_color(c1), new_solid_color(c2));
}

typedef struct noise_texture
{
    texture parent;
    perlin noise;
    float scale;
} noise_texture;

internal float3 noise_texture_value(noise_texture const *this, float2 uv, float3 p)
{
    return 0.5f * (1.0f + sinf(this->scale * p.z + 10 * perlin_turb(&this->noise, p, 7)));
}

__attribute__((overloadable)) internal texture *new_noise_texture(void)
{
    noise_texture *result = malloc(sizeof(noise_texture));
    result->noise = make_perlin();
    result->scale = 0;

    /* setup vtable */
    static struct texture_vtable noise_texture_vtable = { .value = (void*)&noise_texture_value };
    result->parent.vtable = &noise_texture_vtable;

    return (void*)result;
}

__attribute__((overloadable)) internal texture *new_noise_texture(float sc)
{
    noise_texture *result = (void*)new_noise_texture();
    result->scale = sc;
    return (void*)result;
}

aabb.h

#pragma once

typedef struct aabb
{
    float3 min;
    float3 max;
} aabb;

internal bool aabb_hit(aabb const *this, ray const *r, float tmin, float tmax)
{
    #define swap(a, b) ({__typeof__(a) temp = a; a = b; b = temp;})
    for(i32 a = 0; a < 3; ++a)
    {   
        float inv_d  = 1.0f / r->dir[a];
        float t0 = (this->min[a] - r->pos[a]) * inv_d;
        float t1 = (this->max[a] - r->pos[a]) * inv_d;
        if(inv_d < 0.0f)
            swap(t0, t1);
            
        tmin = t0 > tmin ? t0 : tmin;
        tmax = t1 < tmax ? t1 : tmax; 
        if(tmax <= tmin)
            return false;   
    }
    return true;
    #undef swap
}

internal aabb surrounding_box(aabb box0, aabb box1)
{
    float3 small = {fminf(box0.min.x, box1.min.x),
                    fminf(box0.min.y, box1.min.y),
                    fminf(box0.min.z, box1.min.z)};

    float3 big = {fmaxf(box0.max.x, box1.max.x),
                  fmaxf(box0.max.y, box1.max.y),
                  fmaxf(box0.max.z, box1.max.z)};

    return (aabb){small, big};
}

hittable.h

#pragma once

typedef struct material material;
typedef struct hit_record
{
    float3 p;
    float3 normal;
    material *mat_ptr;
    float t;
    float2 uv;
    bool front_face;
} hit_record;

internal inline void hit_record_set_face_normal(hit_record *this, ray const *r, float3 outward_normal)
{
    this->front_face = dot(r->dir, outward_normal) < 0;
    this->normal = this->front_face ? outward_normal : -outward_normal;
}

struct hittable;

struct hittable_vtable
{
    bool(*hit)(struct hittable const *this, ray const *r, float t_min, float t_max, hit_record* rec);
    bool(*bounding_box)(struct hittable const *this, float t0, float t1, aabb *output_box);
};

/* an abstract class */
typedef struct hittable
{
    struct hittable_vtable *vtable; 
} hittable;

typedef struct translate
{
    hittable parent;
    hittable *ptr;
    float3 offset;
} translate;

internal bool translate_hit(translate const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    ray moved_r = {r->pos - this->offset, r->dir, r->time};
    if(!this->ptr->vtable->hit(this->ptr, &moved_r, t_min, t_max, rec))
        return false;

    rec->p += this->offset;
    hit_record_set_face_normal(rec, &moved_r, rec->normal);

    return true;
}

internal bool translate_bounding_box(translate const *this, float t0, float t1, aabb *output_box)
{
    if(!this->ptr->vtable->bounding_box(this->ptr, t0, t1, output_box))
        return false;
    *output_box = (aabb){output_box->min + this->offset,
                         output_box->max + this->offset};
    return true;
}

internal translate make_translate(hittable *ptr, float3 displacement)
{
    translate result = {
        .ptr = ptr,
        .offset = displacement
    };
    /* setup vtable */
    static struct hittable_vtable translate_vtable = {
        .hit = (void*)&translate_hit,
        .bounding_box = (void*)&translate_bounding_box
    };
    result.parent.vtable = &translate_vtable;

    return result;
}

typedef struct rotate_y
{
    hittable parent;
    hittable *ptr;
    aabb bbox;
    float sin_theta;
    float cos_theta;
    bool hasbox;
} rotate_y;

internal bool rotate_y_bounding_box(rotate_y const *this, float t0, float t1, aabb *output_box)
{
    *output_box = this->bbox;
    return this->hasbox;
}

internal bool rotate_y_hit(rotate_y const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    float3 origin = r->pos;
    float3 direction = r->dir;

    origin[0] = this->cos_theta * r->pos[0] - this->sin_theta * r->pos[2];
    origin[2] = this->sin_theta * r->pos[0] + this->cos_theta * r->pos[2];

    direction[0] = this->cos_theta * r->dir[0] - this->sin_theta * r->dir[2];
    direction[2] = this->sin_theta * r->dir[0] + this->cos_theta * r->dir[2];

    ray rotated_ray = {origin, direction, r->time};

    if(!this->ptr->vtable->hit(this->ptr, &rotated_ray, t_min, t_max, rec))
        return false;
        
    float3 p = rec->p;
    float3 normal = rec->normal;

    p[0] =  this->cos_theta * rec->p[0] + this->sin_theta * rec->p[2];
    p[2] = -this->sin_theta * rec->p[0] + this->cos_theta * rec->p[2];


    normal[0] =  this->cos_theta * rec->normal[0] + this->sin_theta * rec->normal[2];
    normal[2] = -this->sin_theta * rec->normal[0] + this->cos_theta * rec->normal[2];

    rec->p = p;
    rec->normal = normal;

    return true;
}

internal rotate_y make_rotate_y(hittable *p, float angle)
{
    float radians = degrees_to_radians(angle);
    rotate_y result = {
        .ptr = p,
        .sin_theta = sinf(radians),
        .cos_theta = cosf(radians),
        .hasbox = p->vtable->bounding_box(p, 0, 1, &result.bbox)
    };

    float3 min = INFINITY;
    float3 max = -INFINITY;

    for(i32 i = 0; i < 2; ++i)
        for(i32 j = 0; j < 2; ++j)
            for(i32 k = 0; k < 2; ++k)
            {
                float x = i * result.bbox.max.x + (1 - i) * result.bbox.min.x;
                float y = i * result.bbox.max.y + (1 - i) * result.bbox.min.y;
                float z = i * result.bbox.max.z + (1 - i) * result.bbox.min.z;

                float newx =  result.cos_theta * x + result.sin_theta * z;
                float newz = -result.sin_theta * x + result.cos_theta * z;

                float3 tester = {newx, y, newz};

                for(i32 c = 0; c < 3; ++c)
                {
                    min[c] = fminf(min[c], tester[c]);
                    max[c] = fmaxf(max[c], tester[c]);
                }
            }

    result.bbox = (aabb){min, max};


    /* setup vtable */
    static struct hittable_vtable rotate_y_vtable = {
        .hit = (void*)&rotate_y_hit,
        .bounding_box = (void*)&rotate_y_bounding_box
    };
    result.parent.vtable = &rotate_y_vtable;
    
    return result;
}

material.h

#pragma once
struct material;

struct material_vtable
{
    bool(*scatter)(struct material const *this, ray const *r_in, hit_record const *rec, float3 *attenuation, ray *scattered);
    float3(*emitted)(struct material const *this, float2 uv, float3 p);
};

typedef struct material
{
    struct material_vtable *vtable;
} material;

internal float3 material_emitted(material const *this, float2 uv, float3 p)
{
    (void)this, (void)uv, (void)p;
    return 0;
}

internal bool material_scatter(material const *this, ray const *r_in, hit_record const *rec, float3 *attenuation, ray *scattered)
{
    return false;
}

typedef struct lambertian
{
    material parent;
    texture *albedo;
} lambertian;

internal bool lambertian_scatter(lambertian const *this, ray const *r_in, hit_record const *rec, float3 *attenuation, ray *scattered)
{
    float3 scatter_direction = rec->normal + random_in_unit_sphere();
    *scattered = (ray){rec->p, scatter_direction, r_in->time};
    *attenuation = this->albedo->vtable->value(this->albedo, rec->uv, rec->p);
    return true;
}

__attribute__((overloadable)) internal material *new_lambertian(texture *color)
{
    lambertian *result = malloc(sizeof(lambertian));
    result->albedo = color;
    
    static struct material_vtable lambertian_vtable = {
        .scatter = (void*)&lambertian_scatter,
        .emitted = &material_emitted
    };
    result->parent.vtable = &lambertian_vtable;
    return (void*)result;
}

__attribute__((overloadable)) internal material *new_lambertian(float3 color)
{
    return new_lambertian(new_solid_color(color));
}

typedef struct metal
{
    material parent;
    float3 albedo;
    float fuzz;
} metal;

internal bool metal_scatter(metal const *this, ray const *r_in, hit_record const *rec, float3 *attenuation, ray *scattered)
{
    float3 reflected = reflect(normalize(r_in->dir), rec->normal);
    *scattered = (ray){rec->p, reflected + this->fuzz * random_in_unit_sphere()};
    *attenuation = this->albedo;
    return (dot(scattered->dir, rec->normal) > 0);
}

internal material *new_metal(float3 color, float fuzz)
{
    metal *result = malloc(sizeof(metal));
    result->albedo = color;
    result->fuzz = fuzz < 1 ? fuzz : 1;

    static struct material_vtable metal_vtable = {
        .scatter = (void*)&metal_scatter,
        .emitted = (void*)&material_emitted
    };
    result->parent.vtable = &metal_vtable;
    return (void*)result;
}

typedef struct dielectric
{
    material parent;
    float ref_idx;
} dielectric;

internal inline float schlick(float cosine, float ref_idx)
{
    float r0 = (1 - ref_idx) / (1 + ref_idx);
    r0 *= r0;
    return r0 + (1 - r0) * powf(1 - cosine, 5);
}

internal bool dielectric_scatter(dielectric const *this, ray const *r_in, hit_record const *rec, float3 *attenuation, ray *scattered)
{
    *attenuation = 1.0f;
    float etai_over_etat = rec->front_face ? 1.0f / this->ref_idx : this->ref_idx;

    float3 unit_direction = normalize(r_in->dir);

    float cos_theta = fminf(dot(-unit_direction, rec->normal), 1.0f);
    float sin_theta = sqrtf(1.0f - cos_theta * cos_theta);
    if(etai_over_etat * sin_theta > 1.0f)
    {
        float3 reflected = reflect(unit_direction, rec->normal);
        *scattered = (ray){rec->p, reflected};
        return true;
    }
    
    float reflect_prob  = schlick(cos_theta, etai_over_etat);
    if(randomf() < reflect_prob)
    {
        float3 reflected = reflect(unit_direction, rec->normal);
        *scattered = (ray){rec->p, reflected};
        return true;
                
    }
    
    float3 refracted = refract(unit_direction, rec->normal, etai_over_etat);
    *scattered = (ray){rec->p, refracted};
    return true;
}

internal material *new_dielectric(float ref_idx)
{
    dielectric *result = malloc(sizeof(dielectric));
    result->ref_idx = ref_idx;

    static struct material_vtable dielectric_vtable = {
        .scatter = (void*)&dielectric_scatter,
        .emitted = &material_emitted
    };
    result->parent.vtable = &dielectric_vtable;
    return (void*)result;
}

typedef struct diffuse_light
{
    material parent;
    texture *emit;
} diffuse_light;

internal float3 diffuse_light_emitted(diffuse_light const *this, float2 uv, float3 p)
{
    return this->emit->vtable->value(this->emit, uv, p);
}

__attribute__((overloadable)) internal material *new_diffuse_light(texture *a)
{
    diffuse_light *result = malloc(sizeof(diffuse_light));
    result->emit = a;

    /* setup vtable */
    static struct material_vtable diffuse_light_vtable = {
        .scatter = &material_scatter,
        .emitted = (void*)&diffuse_light_emitted
    };
    result->parent.vtable = &diffuse_light_vtable;

    return (void*)result;
}

__attribute__((overloadable)) internal material *new_diffuse_light(float3 c)
{
    return new_diffuse_light(new_solid_color(c));
}

typedef struct isotropic
{
    material parent;
    texture *albedo;
} isotropic;

internal bool isotropic_scatter(isotropic const *this, ray const *r, hit_record const *rec, float3 *attenuation, ray *scattered)
{
    *scattered = (ray){rec->p, random_in_unit_sphere(), r->time};
    *attenuation = this->albedo->vtable->value(this->albedo, rec->uv, rec->p);
    return true;
}

__attribute__((overloadable)) internal material *new_isotropic(texture *a)
{
    isotropic *result = malloc(sizeof(isotropic));
    result->albedo = a;

    /* setup vtable */
    static struct material_vtable isotropic_vtable = {
        .scatter = (void*)&isotropic_scatter,
        .emitted = &material_emitted
    };
    result->parent.vtable = &isotropic_vtable;

    return (void*)result;
}

hittable_list.h

#pragma once
typedef struct hittable_list
{
    hittable parent;
    struct
    {
        size_t size;
        size_t capacity;
        hittable **data;
    } objects;
} hittable_list;

internal void hittable_list_reserve(void *self, size_t size)
{
    hittable_list *this = self;

    this->objects.capacity = size;
    this->objects.data = realloc(this->objects.data, size * sizeof(hittable*));
}

internal void hittable_list_add(void *self, void *object)
{
    hittable_list *this = self;
    
    /* if we have not already allocate memory */
    if(NULL == this->objects.data)
    {
        this->objects.data = malloc(sizeof(hittable*));
        this->objects.data[0] = object;
        this->objects.size = 1;
        this->objects.capacity = 1;
    }
    else
    {
        if(++this->objects.size > this->objects.capacity)
        {
            this->objects.capacity = this->objects.size * 2;
            this->objects.data = realloc(this->objects.data, sizeof(hittable*) * this->objects.capacity);
        }

        this->objects.data[this->objects.size - 1] = object;
    }
}

internal bool hittable_list_hit(hittable_list const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    hit_record temp_rec;
    bool hit_anything = false;
    float closest_so_far = t_max;

    for(size_t i = 0; i < this->objects.size; ++i)
    {
        if(this->objects.data[i]->vtable->hit(this->objects.data[i], r, t_min, closest_so_far, &temp_rec))
        {
            hit_anything = true;
            closest_so_far = temp_rec.t;
            *rec = temp_rec;
        }
    }

    return hit_anything;
}

bool hittable_list_bounding_box(hittable_list const *this, float t0, float t1, aabb *output_box)
{
    if(this->objects.size == 0) return false;

    aabb temp_box;
    bool first_box = true;

    for(size_t i = 0; i < this->objects.size; ++i)
    {
        if(!this->objects.data[i]->vtable->bounding_box(this->objects.data[i], t0, t1, &temp_box)) return false;
        *output_box = first_box ? temp_box : surrounding_box(*output_box, temp_box);
        first_box = false;      
    }

    return false;
}


__attribute__((overloadable)) internal hittable_list make_hittable_list(void)
{
    static struct hittable_vtable hittable_list_vtable = { 
        .hit = (void*)&hittable_list_hit,
        .bounding_box = (void*)&hittable_list_bounding_box
    };
    return (hittable_list){ .parent.vtable = &hittable_list_vtable };
}

__attribute__((overloadable)) internal hittable_list make_hittable_list(void *object)
{
    hittable_list result = make_hittable_list();
    hittable_list_add(&result, object);
    return result;
}

sphere.h

#pragma once
typedef struct sphere
{
    hittable parent;
    float3 center;
    float radius;
    material *mat_ptr;  
} sphere;

internal void get_sphere_uv(float3 p, float2 *uv)
{
    float phi = atan2f(p.z, p.z);
    float theta = asinf(p.y);
    (*uv)[0] = 1 - (phi + PI) / (2 * PI);
    (*uv)[1] = (theta + PI / 2) / PI;
}

internal bool sphere_hit(sphere const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    float3 oc = r->pos - this->center;
    float a = length_sqr(r->dir);
    float half_b = dot(oc, r->dir);
    float c = length_sqr(oc) - this->radius * this->radius;
    float discriminant = half_b * half_b - a * c;

    if(discriminant > 0.0f)
    {
        float root = sqrtf(discriminant);
        float temp = (-half_b - root) / a;
        if(temp < t_max && temp > t_min)
        {
            rec->t = temp;
            rec->p = ray_at(r, rec->t);
            float3 outward_normal = (rec->p - this->center) / this->radius;
            hit_record_set_face_normal(rec, r, outward_normal);
            get_sphere_uv((rec->p - this->center) / this->radius, &rec->uv);            
            rec->mat_ptr = this->mat_ptr;
            return true;
        }

        temp = (-half_b + root) / a;
        if(temp < t_max && temp > t_min)
        {
            rec->t = temp;
            rec->p = ray_at(r, rec->t);
            float3 outward_normal = (rec->p - this->center) / this->radius;
            hit_record_set_face_normal(rec, r, outward_normal);
            get_sphere_uv((rec->p - this->center) / this->radius, &rec->uv);
            rec->mat_ptr = this->mat_ptr;
            return true;
        }
    }

    return false;
}

internal bool sphere_bounding_box(sphere const *this, float t0, float t1, aabb *output_box)
{
    *output_box = (aabb){this->center - this->radius, this->center + this->radius};
    return true;
}

sphere make_sphere(float3 center, float radius, material *m)
{
    /* setup vtable */
    static struct hittable_vtable sphere_vtable = {
        .hit = (void*)&sphere_hit,
        .bounding_box = (void*)&sphere_bounding_box
    };

    return (sphere) {
        .center = center,
        .radius = radius,
        .mat_ptr = m,
        .parent.vtable = &sphere_vtable
    };
}

aarect.h

#pragma once

typedef struct xy_rect
{
    hittable parent;
    material *mat_ptr;
    float x0, x1, y0, y1, k;
} xy_rect;

internal bool xy_rect_bounding_box(xy_rect const *this, float t0, float t1, aabb *output_box)
{
    *output_box = (aabb){(float3){this->x0, this->y0, this->k - 0.0001f}, (float3){this->x1, this->y1, this->k - 0.0001f}};
    return true;    
}

internal bool xy_rect_hit(xy_rect const *this, ray const *r, float t0, float t1, hit_record *rec)
{
    float t = (this->k - r->pos.z) / r->dir.z;
    if(t < t0 || t > t1)
        return false;

    float x = r->pos.x + t * r->dir.x;
    float y = r->pos.y + t * r->dir.y;
    if(x < this->x0 || x > this->x1 || y < this->y0 || y > this->y1)
        return false;

    rec->uv = (float2){(x - this->x0) / (this->x1 - this->x0),
                       (y - this->y0) / (this->y1 - this->y0)};
    rec->t = t;
    float3 outward_normal = {0, 0, 1};
    hit_record_set_face_normal(rec, r, outward_normal);
    rec->mat_ptr = this->mat_ptr;
    rec->p = ray_at(r, t);
    return true;
}

internal hittable *new_xy_rect(float x0, float x1, float y0, float y1, float k, material *mat)
{
    xy_rect *result = malloc(sizeof(xy_rect));
    result->x0 = x0;
    result->x1 = x1;
    result->y0 = y0;
    result->y1 = y1;
    result->k = k;
    result->mat_ptr = mat;

    /* setup vtable */
    static struct hittable_vtable xy_rect_vtable = {
        .hit = (void*)&xy_rect_hit,
        .bounding_box = (void*)&xy_rect_bounding_box
    };
    result->parent.vtable = &xy_rect_vtable;
    return (void*)result;
}

typedef struct xz_rect
{
    hittable parent;
    material *mat_ptr;
    float x0, x1, z0, z1, k;
} xz_rect;

internal bool xz_rect_bounding_box(xz_rect const *this, float t0, float t1, aabb *output_box)
{
    *output_box = (aabb){(float3){this->x0, this->k - 0.0001f, this->z0}, (float3){this->x1, this->k - 0.0001f, this->z1}};
    return true;
}

internal bool xz_rect_hit(xz_rect const *this, ray const *r, float t0, float t1, hit_record *rec)
{
    float t = (this->k - r->pos.y) / r->dir.y;
    if(t < t0 || t > t1)
        return false;

    float x = r->pos.x + t * r->dir.x;
    float z = r->pos.z + t * r->dir.z;
    if(x < this->x0 || x > this->x1 || z < this->z0 || z > this->z1)
        return false;

    rec->uv = (float2){(x - this->x0) / (this->x1 - this->x0),
                       (z - this->z0) / (this->z1 - this->z0)};
    rec->t = t;
    float3 outward_normal = {0, 1, 0};
    hit_record_set_face_normal(rec, r, outward_normal);
    rec->mat_ptr = this->mat_ptr;
    rec->p = ray_at(r, t);
    return true;
}

internal hittable *new_xz_rect(float x0, float x1, float z0, float z1, float k, material *mat)
{
    xz_rect *result = malloc(sizeof(xy_rect));
    result->x0 = x0;
    result->x1 = x1;
    result->z0 = z0;
    result->z1 = z1;
    result->k = k;
    result->mat_ptr = mat;

    /* setup vtable */
    static struct hittable_vtable xz_rect_vtable = {
        .hit = (void *) &xz_rect_hit,
        .bounding_box = (void *) &xz_rect_bounding_box
    };
    result->parent.vtable = &xz_rect_vtable;
    return (void*)result;
}

typedef struct yz_rect
{
    hittable parent;
    material *mat_ptr;
    float y0, y1, z0, z1, k;
} yz_rect;

internal bool yz_rect_bounding_box(yz_rect const *this, float t0, float t1, aabb *output_box)
{
    *output_box = (aabb){(float3){this->k - 0.0001f, this->y0, this->z0}, (float3){this->k - 0.0001f, this->y1, this->z1}};
    return true;
}

internal bool yz_rect_hit(yz_rect const *this, ray const *r, float t0, float t1, hit_record *rec)
{
    float t = (this->k - r->pos.x) / r->dir.x;
    if(t < t0 || t > t1)
        return false;

    float y = r->pos.y + t * r->dir.y;
    float z = r->pos.z + t * r->dir.z;
    if(y < this->y0 || y > this->y1 || z < this->z0 || z > this->z1)
        return false;

    rec->uv = (float2){(y - this->y0) / (this->y1 - this->y0),
                       (z - this->z0) / (this->z1 - this->z0)};
    rec->t = t;
    float3 outward_normal = {1, 0, 0};
    hit_record_set_face_normal(rec, r, outward_normal);
    rec->mat_ptr = this->mat_ptr;
    rec->p = ray_at(r, t);
    return true;
}

internal hittable *new_yz_rect(float y0, float y1, float z0, float z1, float k, material *mat)
{
    yz_rect *result = malloc(sizeof(yz_rect));
    result->y0 = y0;
    result->y1 = y1;
    result->z0 = z0;
    result->z1 = z1;
    result->k = k;
    result->mat_ptr = mat;

    /* setup vtable */
    static struct hittable_vtable yz_rect_vtable = {
        .hit = (void*)&yz_rect_hit,
        .bounding_box = (void*)&yz_rect_bounding_box,
    };
    result->parent.vtable = &yz_rect_vtable;
    return (void*)result;
}

box.h

#pragma once
typedef struct box
{
    hittable parent;
    float3 box_min;
    float3 box_max;
    hittable_list sides;
} box;

internal bool box_hit(box const *this, ray const *r, float t0, float t1, hit_record *rec)
{
    return hittable_list_hit(&this->sides, r, t0, t1, rec);
}

internal bool box_bounding_box(box const *this, float t0, float t1, aabb *output_box)
{
    *output_box = (aabb){this->box_min, this->box_max};
    return true;
}

internal box make_box(float3 p0, float3 p1, material *mat)
{
    box result = {
        .box_min = p0,
        .box_max = p1,
        .sides = make_hittable_list()
    };

    /* add sides to the result box */
    hittable_list_add(&result.sides, new_xy_rect(p0.x, p1.x, p0.y, p1.y, p1.z, mat));
    hittable_list_add(&result.sides, new_xy_rect(p0.x, p1.x, p0.y, p1.y, p0.z, mat));

    hittable_list_add(&result.sides, new_xz_rect(p0.x, p1.x, p0.z, p1.z, p1.y, mat));
    hittable_list_add(&result.sides, new_xz_rect(p0.x, p1.x, p0.z, p1.z, p0.y, mat));
    
    hittable_list_add(&result.sides, new_yz_rect(p0.y, p1.y, p0.z, p1.z, p1.x, mat));
    hittable_list_add(&result.sides, new_yz_rect(p0.y, p1.y, p0.z, p1.z, p0.x, mat));

    /* setup vtable */
    static struct hittable_vtable box_vtable = {
        .hit = (void*)box_hit,
        .bounding_box = (void*)box_bounding_box
    };
    result.parent.vtable = &box_vtable;
    
    return result;
}

constant_medium.h

#pragma once
typedef struct constant_medium
{
    hittable parent;
    hittable *boundary;
    material *phase_function;
    float neg_inv_density;
} constant_medium;

internal bool constant_medium_bounding_box(constant_medium const *this, float t0, float t1, aabb *output_box)
{
    return this->boundary->vtable->bounding_box(this->boundary, t0, t1, output_box);
}

internal bool constant_medium_hit(constant_medium const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    hit_record rec1, rec2;

    if(!this->boundary->vtable->hit(this->boundary, r, -INFINITY, INFINITY, &rec1))
        return false;

    if(!this->boundary->vtable->hit(this->boundary, r, rec1.t + 0.0001f, INFINITY, &rec2))
        return false;

    if(rec1.t < t_min) rec1.t = t_min;
    if(rec2.t > t_max) rec2.t = t_max;

    if(rec1.t >= rec2.t)
        return false;

    if(rec1.t < 0)
        rec1.t = 0;

    float ray_length = length(r->dir);
    float distance_inside_boundary = (rec2.t - rec1.t) * ray_length;
    float hit_distance = this->neg_inv_density * logf(randomf());

    if(hit_distance > distance_inside_boundary)
        return false;

    rec->t = rec1.t + hit_distance /  ray_length;
    rec->p = ray_at(r, rec->t);
    rec->normal = (float3){1, 0, 0};
    rec->front_face = false;
    rec->mat_ptr = this->phase_function;

    return true;
}

__attribute__((overloadable)) internal constant_medium make_constant_medium(hittable *bounds, float density, texture *a)
{
    constant_medium result = {
        .boundary = bounds,
        .neg_inv_density = -1 / density,
        .phase_function = new_isotropic(a)
    };

    /* setup vtable */
    static struct hittable_vtable constant_medium_vtable = {
        .hit = (void*)&constant_medium_hit,
        .bounding_box = (void*)&constant_medium_bounding_box
    };
    result.parent.vtable = &constant_medium_vtable;

    return result;
}

__attribute__((overloadable)) internal constant_medium make_constant_medium(hittable *bounds, float density, float3 color)
{
    return make_constant_medium(bounds, density, new_solid_color(color));
}

bvh.h

#pragma once

typedef struct bvh_node
{
    hittable parent;
    hittable *left;
    hittable *right;
    aabb box;
} bvh_node;

bool bvh_node_bounding_box(bvh_node const *this, float t0, float t1, aabb *output_box)
{
    *output_box = this->box;
    return true;
}

internal bool bvh_node_hit(bvh_node const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    if(!aabb_hit(&this->box, r, t_min, t_max))
        return false;

    bool hit_left = this->left->vtable->hit(this->left, r, t_min, t_max, rec);
    bool hit_right = this->right->vtable->hit(this->right, r, t_min, hit_left ? rec->t : t_max, rec);

    return hit_left || hit_right;   
}

internal inline int box_compare(hittable ** a, hittable ** b, i32 axis)
{
    aabb box_a;
    aabb box_b;

    if(!(*a)->vtable->bounding_box(*a, 0, 0, &box_a)
    || !(*b)->vtable->bounding_box(*b, 0, 0, &box_b))
    {
            fprintf(stderr, "No bounding box in bvh_node constructor.\n");  
    }

    return box_a.min[axis] == box_b.min[axis] ? 0 : box_a.min[axis] < box_b.min[axis] ? 1 : -1;
}

internal int box_x_compare(hittable ** a, hittable ** b)
{
    return box_compare(a, b, 0);
}

internal int box_y_compare(hittable ** a, hittable ** b)
{
    return box_compare(a, b, 1);
}

internal int box_z_compare(hittable ** a, hittable ** b)
{
    return box_compare(a, b, 2);
}

__attribute__((overloadable)) internal hittable *new_bvh_node(hittable **objects, size_t start, size_t end, float time0, float time1)
{
    /* allocate memory for the result */
    bvh_node *result = calloc(sizeof(bvh_node), 1);
    
    i32 axis = (i32)randomf(0,3);
    /* stupid qsort */ int(*comparator)(hittable **, hittable **) = axis == 0 ? box_x_compare
                                                                                            : axis == 1 ? box_y_compare
                                                                                            : box_z_compare;
    size_t object_span = end - start;
    if(object_span == 1)
    {
        result->left = result->right = objects[start];
    }
    else if(object_span == 2)
    {
        if(comparator(&objects[start], &objects[start + 1]) == 1)
        {
            result->left = objects[start];
            result->right = objects[start + 1];
        }
        else
        {
            result->left = objects[start + 1];
            result->right = objects[start];
        }
    }
    else
    {
        qsort(objects + start, object_span, sizeof(hittable*), (void*)comparator);

        size_t mid = start + object_span / 2;
        result->left = new_bvh_node(objects, start, mid, time0, time1);
        result->right = new_bvh_node(objects, mid, end, time0, time1);
    }

    aabb box_left = {0}, box_right = {0};

    if(!result->left->vtable->bounding_box(result->left, time0, time1, &box_left)
    || !result->right->vtable->bounding_box(result->right, time0, time1, &box_right))
    {
        fprintf(stderr, "No bounding box in bvh_node constructor.\n");  
    }
    result->box = surrounding_box(box_left, box_right);

    static struct hittable_vtable bvh_node_vtable = {
        .hit = (void*)&bvh_node_hit,
        .bounding_box = (void*)&bvh_node_bounding_box
    };
    result->parent.vtable = &bvh_node_vtable;
    
    return (void*)result;
}

__attribute__((overloadable)) internal hittable *new_bvh_node(hittable *list_in, float time0, float time1)
{
    hittable_list *list = (void*)list_in;
    return new_bvh_node(list->objects.data, 0, list->objects.size, time0, time1);   
}

random.cc

#include <random>
#include <ctime>

/* NOTE: this is in c++ because rand is not very good */
extern "C" float random_float()
{
    static thread_local std::mt19937 engine {};
    std::uniform_real_distribution<float> dist{0.0f, 1.0f};
    return dist(engine);    
}

moving_sphere.h

#pragma once
typedef struct moving_sphere
{
    hittable parent;
    float3 center0, center1;
    float time0, time1; 
    float radius;
    material *mat_ptr;  
} moving_sphere;

internal float3 moving_sphere_center(moving_sphere const *this, float time)
{
    return this->center0 + ((time - this->time0) / (this->time1 - this->time0)) * (this->center1 - this->center0);
}

internal bool moving_sphere_hit(moving_sphere const *this, ray const *r, float t_min, float t_max, hit_record *rec)
{
    float3 oc = r->pos - moving_sphere_center(this, r->time);
    float a = length_sqr(r->dir);
    float half_b = dot(oc, r->dir);
    float c = length_sqr(oc) - this->radius * this->radius;
    float discriminant = half_b * half_b - a * c;
    
    if(discriminant > 0.0f)
    {
        float root = sqrtf(discriminant);
        float temp = (-half_b - root) / a;
        if(temp < t_max && temp > t_min)
        {
            rec->t = temp;
            rec->p = ray_at(r, rec->t);
            float3 outward_normal = (rec->p - moving_sphere_center(this, r->time)) / this->radius;
            hit_record_set_face_normal(rec, r, outward_normal);
            rec->mat_ptr = this->mat_ptr;
            return true;
        }

        temp = (-half_b + root) / a;
        if(temp < t_max && temp > t_min)
        {
            rec->t = temp;
            rec->p = ray_at(r, rec->t);
            float3 outward_normal = (rec->p - moving_sphere_center(this, r->time)) / this->radius;
            hit_record_set_face_normal(rec, r, outward_normal);
            rec->mat_ptr = this->mat_ptr;
            return true;
        }
    }

    return false;
}

bool moving_sphere_bounding_box(moving_sphere const *this, float t0, float t1, aabb *output_box)
{
    aabb box0 = {moving_sphere_center(this, t0) - this->radius, moving_sphere_center(this, t0) + this->radius};
    aabb box1 = {moving_sphere_center(this, t1) - this->radius, moving_sphere_center(this, t1) + this->radius};
    *output_box = surrounding_box(box0, box1);
    return true;
}

hittable *new_moving_sphere(float3 center0, float3 center1, float t0, float t1, float radius, void *m)
{
    moving_sphere *result = malloc(sizeof(moving_sphere));

    /* assign the function paremeters to their respective struct members */
    result->center0 = center0;
    result->center1 = center1;
    result->time0 = t0;
    result->time1 = t1;
    result->radius = radius;
    result->mat_ptr = m;
    
    /* setup vtable */
    static struct hittable_vtable moving_sphere_vtable = {
        .hit = (void*)&moving_sphere_hit,
        .bounding_box  = (void*)&moving_sphere_bounding_box
    };
    result->parent.vtable = &moving_sphere_vtable;

    return (void*)result;
}

here is what the program generates: enter image description here

\$\endgroup\$
3
  • 1
    \$\begingroup\$ The file moving_sphere.h is missing. To be honest, this doesn't look completely like C. More like a cross between the 2, my C compiler really doesn't like some of this. \$\endgroup\$
    – pacmaninbw
    Sep 22, 2020 at 1:14
  • \$\begingroup\$ this is c it just uses a lot of clang and gcc extensions gcc.gnu.org/onlinedocs/gcc/C-Extensions.html#C-Extensions clang.llvm.org/docs/… \$\endgroup\$
    – nullptr
    Sep 22, 2020 at 1:33
  • \$\begingroup\$ You mention that you are using gcc. I also use gcc, where my installation includes its own types.h, but here you are including an alternate types.h, where there are some interesting, but otherwise confusing definitions. eg.: #define internal static. Why not just use static? And I am not sure what you are gaining by not just using the default types.h? \$\endgroup\$
    – ryyker
    Oct 8, 2020 at 13:03

1 Answer 1

4
\$\begingroup\$

Choice of programming language

It's nice to learn and experience multiple programming languages. It can even be a nice challenge to try to program something in a simpler language like C. However, what you are doing is mixing C and C++ in one program, and even the C code is not standard C, it's using GNU extensions. Mixing things like that is something I would not recommend for serious projects. Some of the GNU extensions you use are just enabling features in C in a non-portable way, that you could get in C++ in a standard, portable way, like function overloading. So either I suggest you go all-in on your challenge to code in C, and do it in portable, idiomatic C11, or just write everything in C++.

Avoid typedef'ing standard types

I would recommend against making typedefs for standard types like uint8_t. Yes, it might save typing a few characters here and there, but now you have to worry about including the header file that defines them. And while it make be fine if you wrote all the code yourself, imagine that you are writing a library. You don't want to declare types in the global namespace, as this might conflict with typedefs from the main application or other libraries.

Another problem is that some of the typedefs might be outright lies on some platforms: neither C nor C++ guarantee that a float is 32 bits and a double is 64 bits.

Related to this, why did you #define internal static? Why not use static directly in the code?

Consider using a vector math library

Trying to implement things from scratch is a nice way to learn about how computers work. However, I recommend you don't try to write everything from scratch, only focus on one subject, and use (standard) libraries for the rest where possible. You are already using SDL and the C standard library. I also recommend you look for a vector mathematics library to use. This frees you from implementing a lot of the tedious maths programming, and allows you to focus more on the ray tracing algorithm itself. If you were writing in C++, I strongly recommend you have a look at GLM. For C a possible choice is the GNU Scientific Library.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.