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 vertex {
    vec4 position;
    vec4 normal;
};
struct triangle {
    vertex v1;
    vertex v2;
    vertex v3;
    aabb bb;
};

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 next_id;
    int primitive_id;
    uint[2] padding;
};

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

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

    uint spheres_count;
    // uint triangles_count;
    bool enable_dof;
    bool debug_bvh;
    uint padding;
};

layout(set = 0, binding = 1) buffer geometry {
    sphere[] spheres;
    // triangle[] triangles;
};

layout(set = 0, binding = 2) buffer bvh {
    bvh_node[] nodes;
};

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

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

uint pcg_hash(inout uint seed)
{
    uint state = seed;
    seed *= 747796405u + 2891336453u;
    uint word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

float rand(inout uint seed)
{
    return float(pcg_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);
}

vec3 diffuse_reflection_cosine_weighted(inout uint seed) {
    float r = rand(seed);
    float radius  = sqrt(r);
    float theta = 2.0 * PI * rand(seed);

    return vec3(cos(theta) * radius, sin(theta) * radius, sqrt(1 - r));
}

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) {
    const vec3 invdir = 1.0 / r.direction;
    const vec3 f = (b.max.xyz - r.origin.xyz) * invdir;
    const vec3 n = (b.min.xyz - r.origin.xyz) * invdir;

    const vec3 tmax = max(f, n);
    const vec3 tmin = min(f, n);

    const float t1 = min(tmax.x, min(tmax.y, tmax.z));
    const float t0 = max(max(tmin.x, max(tmin.y, tmin.z)), 0.0f);

    return t1 >= t0;
}

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

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

    return hit_count;
}

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_triangle(triangle tri, ray r, float min_t, float max_t, out hit_info info) {
    vec3 v2v1 = tri.v2.position.xyz - tri.v1.position.xyz;
    vec3 v3v1 = tri.v3.position.xyz - tri.v1.position.xyz;

    vec3 originv1 = r.origin - tri.v1.position.xyz;
    vec3 normal = cross(v2v1, v3v1);
    vec3 q = cross(originv1, r.direction);
    float d = 1 / dot(r.direction, normal);
    float u = d * dot(-q, v3v1);
    float v = d * dot(q, v2v1);
    float t = d * dot(-normal, originv1);

    if (u < 0 || u > 1 || v < 0|| u + v > 1) {
        return false;
    }

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

    vec3 outward_normal = normal;
    info.front_face = dot(outward_normal, r.direction) < 0.0;
    info.normal = info.front_face ? outward_normal : -outward_normal;
    info.mat = material(vec4(1.0), vec4(0.0), 1.0, 1.5, LAMBERTIAN_MATERIAL, 0);

    return true;
}

bool hit_node(ray r, float min_t, float max_t, out hit_info info) {
    hit_info temp_info;
    float closest_t = max_t;
    bool hit = false;
    int id = 0;

    while(id != -1) {
        // Leafs
        if (nodes[id].primitive_id != -1) {
            if (hit_sphere(spheres[nodes[id].primitive_id], r, min_t, closest_t, temp_info)) {
                info = temp_info;
                closest_t = temp_info.t;
                hit = true;
            }
            // if (hit_triangle(triangles[nodes[id].primitive_id], r, min_t, closest_t, temp_info)) {
            //     info = temp_info;
            //     closest_t = temp_info.t;
            //     hit = true;
            // }
            id = nodes[id].next_id;
        } else { // Internal nodes
            if (hit_aabb(nodes[id].bb, r, min_t, closest_t)) {
                id++;
            } else {
                id = nodes[id].next_id;
            }
        }
    }

    return hit;
}

bool hit_spheres(ray r, float min_t, float max_t, out hit_info info) {
    return hit_node(r, min_t, max_t, info);
}

bool hit_triangles(ray r, float min_t, float max_t, out hit_info info) {
    return hit_node(r, min_t, max_t, info);
}

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) {
    // while (hit_triangles(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) {
            lambertian_scattering(info.mat, current_ray, seed, info, attenuation, scattered);
            final_color += info.mat.emissive * throughput;
            throughput *= attenuation;

            // float NdotR = dot(info.normal, scattered.direction);
            // float PDF = NdotR / PI;
            // vec4 BRDF = attenuation / PI;
            // throughput = (throughput * NdotR / PDF) * BRDF;
            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;
    // return final_color;
}

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 lens_disk = vec2(0.0);
    if (enable_dof) {
        lens_disk = rand_disk * cam.lens_radius;
    }

    vec3 offset = lens_disk.x * cam.right.xyz + lens_disk.y * cam.up.xyz;
    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.025, 0.025, 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)));
}