layout(local_size_x = 8, local_size_y = 8, local_size_z = 1) in;

struct ray {
    vec3 origin;
    vec3 direction;
};

const float PI = 3.141592653589793;
const uint LAMBERTIAN_MATERIAL  = 0;
const uint METAL_MATERIAL       = 1;
const uint DIELECTRIC_MATERIAL  = 2;
const uint EMISSIVE_MATERIAL    = 3;

const vec3 light_position = vec3(10, -25, 0);

struct material {
    vec4 albedo;
    vec4 emissive;
    float fuzz;
    float ior;
    uint type;
    uint padding;
};

struct hit_info {
    float t;
    vec3 point;
    vec3 normal;
    material mat;
    bool front_face;
};

struct aabb {
    vec4 min;
    vec4 max;
};

struct sphere {
    vec4 position;
    material mat;
    aabb bb;
    float radius;
    float[3] padding;
};

struct camera {
    vec4 position;

    vec4 forward;
    vec4 up;
    vec4 right;

    vec4 horizontal;
    vec4 vertical;
    vec4 first_pixel;

    float lens_radius;
    float[3] padding;
};

struct bvh_node {
    aabb bb;
    int left;
    int sphere_id;
    uint[2] padding;
};

layout(set = 0, binding = 0) buffer Inputs {
    vec4 sky_color;
    vec4 ground_color;
    camera cam;

    sphere[512] spheres;
    bvh_node[1024] nodes;

    uint max_bounce;
    uint width;
    uint height;
    uint sample_index;

    uint spheres_count;
    bool debug_bvh;
    uint[2] padding;
};

layout(set = 0, binding = 1, rgba32f) uniform writeonly image2D output_image;
layout(set = 0, binding = 2, rgba32f) uniform readonly image2D accumulation_image;

vec3 at(ray r, float t) {
    return r.origin + r.direction * t;
}


uint wang_hash(inout uint seed)
{
    seed = (seed ^ 61) ^ (seed >> 16);
    seed *= 9;
    seed = seed ^ (seed >> 4);
    seed *= 0x27d4eb2d;
    seed = seed ^ (seed >> 15);
    return seed;
}

float rand(inout uint seed)
{
    return float(wang_hash(seed)) / 4294967296.0;
}

vec3 unit_vec(inout uint seed) {
    float z  = rand(seed) * 2.0 - 1.0;
    float a = 2.0 * PI * rand(seed);
    float r = sqrt(1.0 - z * z);

    return vec3(cos(a) * r, sin(a) * r, z);
}

vec2 disk_vec(inout uint seed) {
    float z  = rand(seed) * 2.0 - 1.0;
    float a = 2.0 * PI * rand(seed);
    float r = sqrt(1.0 - z * z);

    return vec2(cos(a) * r, sin(a) * r);
}

bool hit_aabb(aabb b, ray r, float min_t, float max_t) {
    for (int axis_index = 0; axis_index < 3; axis_index++) {
        float inv_dir = 1.f / r.direction[axis_index];
        float t0 = (b.min[axis_index] - r.origin[axis_index]) * inv_dir;
        float t1 = (b.max[axis_index] - r.origin[axis_index]) * inv_dir;

        if (inv_dir < 0.f) {
            float temp = t0;
            t0 = t1;
            t1 = temp;
        }

        min_t = t0 > min_t ? t0 : min_t;
        max_t = t1 < max_t ? t1 : max_t;

        if (max_t <= min_t) {
            return false;
        }
    }

    return true;
}

uint hit_aabbs(ray r, float min_t, float max_t) {
    uint hit_count = 0;
    for (uint idx = 0; idx <= spheres_count * 2; idx++) {
        hit_count += uint(hit_aabb(nodes[idx].bb, r, min_t, max_t));

        if (nodes[idx].sphere_id >= 0) {
            hit_count += uint(hit_aabb(spheres[nodes[idx].sphere_id].bb, r, min_t, max_t));
        }
    }

    return hit_count;
}

// bool hit_node(ray r, float min_t, float max_t, out hit_info info) {
//     if (!hit_aabb(nodes[0].bb, r, min_t, max_t)) {
//         return false;
//     }
//     return true;
// 
//     int stack_id = -1;
//     int stack[64];
//     bvh_node current_node = nodes[0];
// 
//     while(true) {
//         bool hit_left = hit_aabb(nodes[current_node.left - 1].bb, r, min_t, max_t);
//         bool hit_right = hit_aabb(nodes[current_node.right - 1].bb, r, min_t, max_t);
// 
//         if (hit_left && hit_right) {
//             stack_id++;
//             stack[stack_id] = current_node.left - 1;
// 
//             current_node = nodes[current_node.right - 1];
//         } else if (hit_left) {
//             current_node = nodes[current_node.left - 1];
//         } else if (hit_right) {
//             current_node = nodes[current_node.right - 1];
//         } else {
//             if (stack_id >= 0) {
//                 current_node = nodes[stack[stack_id]];
//                 stack_id--;
//             } else {
//                 return false;
//             }
//         }
// 
//         // Reached a leaf
//         if (current_node.left == -1 && current_node.right == -1) {
//             return true;
//         }
//     }
// 
// }

bool hit_sphere(sphere s, ray r, float min_t, float max_t, out hit_info info) {
    vec3 origin_center = r.origin - s.position.xyz;

    float a = dot(r.direction, r.direction);
    float half_b = dot(r.direction, origin_center);
    float c = dot(origin_center, origin_center) - s.radius * s.radius;
    float discriminant = half_b * half_b - a * c;

    if (discriminant < 0) {
        return false;
    }

    float sqrt_discriminant = sqrt(discriminant);
    float root = (-half_b - sqrt_discriminant) / a;
    if (root < min_t || root > max_t) {
        root = (-half_b + sqrt_discriminant) / a;

        if (root < min_t || root > max_t) {
            return false;
        }
    }

    info.t = root;
    info.point = at(r, info.t);

    vec3 outward_normal = (info.point - s.position.xyz) / s.radius;
    info.front_face = dot(outward_normal, r.direction) < 0.0;
    info.normal = info.front_face ? outward_normal : -outward_normal;
    info.mat = s.mat;

    return true;
}

bool hit_spheres(ray r, float min_t, float max_t, out hit_info info) {
    float closest_t = max_t;
    hit_info temp_info;
    bool hit = false;
    for (uint sphere_index = 0; sphere_index <= spheres_count; sphere_index++) {
        if (spheres[sphere_index].radius > 0.0 && hit_sphere(spheres[sphere_index], r, min_t, closest_t, temp_info)) {
            info = temp_info;
            closest_t = temp_info.t;
            hit = true;
        }
    }

    return hit;
}

bool lambertian_scattering(material mat, ray incident, inout uint seed, hit_info info, out vec4 attenuation, out ray scattered) {
    attenuation = mat.albedo;
    scattered = ray(info.point, normalize(info.normal + unit_vec(seed)));

    return true;
}

bool metal_scattering(material mat, ray incident, inout uint seed, hit_info info, out vec4 attenuation, out ray scattered) {
    vec3 reflected = reflect(normalize(incident.direction), info.normal);

    attenuation = mat.albedo;
    scattered = ray(info.point, reflected + mat.fuzz * unit_vec(seed)); 

    return (dot(scattered.direction, info.normal) > 0.0);
}

float reflectance(float cosine, float eta) {
    // Use Schlick's approximation for reflectance.
    float r0 = (1 - eta) / (1 + eta);
    r0 = r0 * r0;
    return r0 + (1 - r0) * pow((1 - cosine), 5);
}

bool dielectric_scattering(material mat, ray incident, inout uint seed, hit_info info, out vec4 attenuation, out ray scattered) {
    attenuation = vec4(1.0);

    vec3 unit_dir = normalize(incident.direction);
    float cos_theta = min(dot(-unit_dir, info.normal), 1.0);
    float sin_theta = sqrt(1.0 - cos_theta*cos_theta);

    float eta = info.front_face ? 1.0 / mat.ior : mat.ior;
    vec3 direction;

    if (eta * sin_theta > 1.0 || reflectance(cos_theta, eta) > rand(seed)) {
        direction = reflect(unit_dir, info.normal);
    } else {
        direction = refract(unit_dir, info.normal, eta);
    }

    scattered = ray(info.point, direction); 

    return true;
}

vec4 ray_color(ray r, uint seed) {
    hit_info info;
    uint bounces = 0;

    vec4 final_color = vec4(0.0);
    vec4 throughput = vec4(1.0);
    ray current_ray = r;
    while (hit_spheres(current_ray, 0.001, 1e15, info) && bounces < max_bounce) {
        ray scattered;
        vec4 attenuation = vec4(0.0);
        if ((info.mat.type == LAMBERTIAN_MATERIAL && lambertian_scattering(info.mat, current_ray, seed, info, attenuation, scattered)) ||
        (info.mat.type == METAL_MATERIAL && metal_scattering(info.mat, current_ray, seed, info, attenuation, scattered)) ||
        (info.mat.type == DIELECTRIC_MATERIAL && dielectric_scattering(info.mat, current_ray, seed, info, attenuation, scattered)) ||
        info.mat.type == EMISSIVE_MATERIAL) {
            final_color += info.mat.emissive * throughput;
            throughput *= attenuation;
            current_ray = scattered;

            // Russian Roulette
            // As the throughput gets smaller, the ray is more likely to get terminated early.
            // Survivors have their value boosted to make up for fewer samples being in the average.
            // {
            //     float p = max(throughput.r, max(throughput.g, throughput.b));
            //     if (rand(seed) > p) {
            //         break;
            //     }

            //     // Add the energy we 'lose' by randomly terminating paths
            //     throughput *= 1.0 / p;
            // }
        }
        bounces++;
    }

    return throughput;
}

void main() {
    uint x = gl_GlobalInvocationID.x;
    uint y = gl_GlobalInvocationID.y;

    vec2 uv = vec2(float(x) / (width - 1), 1.0 - (float(y) / (height - 1)));

    uint seed = uint(x * uint(1973) + y * uint(9277) + sample_index * uint(26699)) | uint(1);
    // vec2 rand_disk = disk_vec(seed);
    // vec2 jittered_uvs = uv + vec2(rand(seed) / (width - 1), rand(seed) / (height - 1));
    vec2 jittered_uvs = uv;
    // vec2 lens_disk = rand_disk * cam.lens_radius;
    // vec3 offset = lens_disk.x * cam.right.xyz + lens_disk.y * cam.up.xyz;
    vec3 offset = vec3(0.0);
    vec4 proj_plane_pos = cam.first_pixel + jittered_uvs.x * cam.horizontal + jittered_uvs.y * cam.vertical;

    ray r = ray(cam.position.xyz + offset, proj_plane_pos.xyz - cam.position.xyz - offset);

    vec4 out_color = ray_color(r, seed);

    if (debug_bvh) {
        out_color -= hit_aabbs(r, 0.001, 1e15) * vec4(0.0, 0.05, 0.05, 0.0);
    }

    float spp_scale = (1.0 / sample_index);
    vec4 acc_color = imageLoad(accumulation_image, ivec2(x, y));
    imageStore(output_image, ivec2(x, y), vec4(mix(acc_color, out_color, spp_scale)));
}