#define SAMPLES 8
#define BOUNCES 8

// sphere of size ra centered at point ce
vec2 sphIntersect(in vec3 ro, in vec3 rd, in vec3 ce, float ra)
{
    vec3 oc = ro - ce;
    float b = dot( oc, rd );
    vec3 qc = oc - b*rd;
    float h = ra*ra - dot( qc, qc );
    if( h<0.0 ) return vec2(-1.0); // no intersection
    h = sqrt( h );
    return vec2( -b-h, -b+h );
}

// disk center c, normal n, radius r
float diskIntersect( in vec3 ro, in vec3 rd, vec3 c, vec3 n, float r )
{
    vec3  o = ro - c;
    float t = -dot(n,o)/dot(rd,n);
    vec3  q = o + rd*t;
    return (dot(q,q)<r*r) ? t : -1.0;
}

struct Mat {
    vec3 emission;
    vec3 albedo;
};

const Mat WHITE = Mat(vec3(0), vec3(0.89));
const Mat LIGHT = Mat(vec3(4, 2, 0) * 32.0, vec3(0));
const Mat LIGHT2 = Mat(vec3(0, 4, 2) * 32.0, vec3(0));

struct hit {
    vec3 nor;
    float t;
    Mat mat;
} Hit;

void intersect_sp(in vec3 ro, in vec3 rd, in vec3 o, in float r, in Mat mat)
{
    vec2 hit = sphIntersect(ro, rd, o, r);
    float d = hit.x;

    if (hit.x < 1e-3)
        d = hit.y;

    if (d < Hit.t && d > 1e-3)
    {
        Hit.t = d;
        Hit.nor = normalize(ro + rd * d - o);
        Hit.mat = mat;
    }
}


void intersect_disk(in vec3 ro, in vec3 rd, in vec3 o, in vec3 n, in float r, in Mat mat)
{
    float d = diskIntersect(ro, rd, o, n, r);

    if (d < Hit.t && d > 1e-3)
    {
        Hit.t = d;
        Hit.nor = n;
        Hit.mat = mat;
    }
}


void intersect(in vec3 ro, in vec3 rd)
{
    Hit.t = 1e3;

    intersect_sp(ro, rd, vec3(0,1,0), 1.0, WHITE);
    intersect_sp(ro, rd, vec3(2,1,0), 1.0, WHITE);
    intersect_sp(ro, rd, vec3(2,3,-0.5), 0.5, LIGHT);
    intersect_sp(ro, rd, vec3(-2,4,+0.5), 0.5, LIGHT2);
    intersect_sp(ro, rd, vec3(0,3,0), 1.0, WHITE);
    intersect_disk(ro, rd, vec3(0,0,0), vec3(0,1,0), 4.0, WHITE);
}

vec3 cosine_hemisphere(in float ra)
{
    float u1 = gold_noise();
    float u2 = gold_noise();

    float r = sqrt(u1) * ra;
    float theta = 2.0 * PI * u2;

    float x = r * cos(theta);
    float y = r * sin(theta);

    return vec3(x, y, sqrt(1.0 - u1));
}

// project vector b onto vector a
vec3 project(in vec3 a, in vec3 b) {
    return a * (dot(a, b) / dot(a, a));
}

// construct a 3D coordinate system with the input up being the "upwards" facing vector
// which will be directly stored in w.
// Mathematically this function will create two non linear vectors of up and generate an orthonormal
// basis using gram-schmidt.
// This function assumes "up" being already normalized
void construct_orthonormal_basis(in vec3 up, out vec3 u, out vec3 v, out vec3 w) {
    w = up;

    vec3 n2 = normalize(cross(w, vec3(0, 1, 1)));   // build perpendicular vector from w
    vec3 n3 = cross(w, n2); // create 2nd vector perpendicular to w and n2

    // gram schmidt
    u = n2 - (project(w, n2));
    v = n3 - project(w, n3) - project(u, n3);
}

vec3 gen_diffuse(in vec3 normal, in vec3 incident) {
    vec3 hemisphere = cosine_hemisphere(1.0);

    vec3 u, v, w;
    construct_orthonormal_basis(normal, u, v, w);

    return u * hemisphere.x + v * hemisphere.y + w * hemisphere.z;
}

vec3 gen_reflection(in vec3 normal, in vec3 incident, in float roughness) {
    vec3 hemisphere = cosine_hemisphere(roughness);

    vec3 u, v, w;
    construct_orthonormal_basis(reflect(incident, normal), u, v, w);

    return u * hemisphere.x + v * hemisphere.y + w * hemisphere.z;
}

float schlick(in float cosTheta, in float R0) {
    return R0 + (1.0 - R0) * pow(1.0 - cosTheta, 5.0);
}

float R0(in float ior) {
    float a = (ior - 1.0);
    float b = (ior + 1.0);
    return (a * a) / (b * b);
}

vec3 trace_path(in vec3 ro, in vec3 rd)
{
    vec3 rad = vec3(1);
    vec3 acc = vec3(0);

    for (int i = 0; i < BOUNCES; i++)
    {
        intersect(ro, rd);
        if (Hit.t < 1e3)
        {
            float cosTheta = abs(dot(rd, Hit.nor));
            float pdf = cosTheta * INV_PI;

            acc += rad * Hit.mat.emission * pdf;

            ro = ro + rd * Hit.t;

            if (schlick(cosTheta, R0(1.450)) < gold_noise())
            {
                rd = gen_diffuse(Hit.nor, rd);
                rad *= Hit.mat.albedo;
            } else
            {
                rd = gen_reflection(Hit.nor, rd, 0.05);
            }
        } else {
            acc += rad * pow(texture(iChannel1, rd).rgb, vec3(2.0)) * 8.0;
            break;
        }
    }

    return acc;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    init_random_state(fragCoord.xy, iTime);

    vec2 taa_off = vec2(gold_noise(), gold_noise()) * 1.5 - 0.5;
    vec2 uv = (fragCoord + taa_off - 0.5 * iResolution.xy) / iResolution.y;

    vec3 col = vec3(0);

    float s = sin(iTime);
    float c = cos(iTime);
    mat2 r = mat2(c, -s, s, c);

    vec3 rd = normalize(vec3(uv, 1.0));
    vec3 ro = vec3(0,2,-8.0);
    float t = 0.01;

    rd.xz *= r;
    ro.xz *= r;

    for (int i = 0; i < SAMPLES; i++) {
        col += trace_path(ro, rd);
    }
    col /= float(SAMPLES);

    vec3 prev = texture(iChannel0, fragCoord/iResolution.xy).rgb;

    fragColor = vec4(mix(prev, col, 0.3) ,1.0);
}