#version 300 es

precision highp float;

out vec4 outColor;

// --- inversive geometry ---

struct vecInv {
    vec3 sp;
    vec2 lt;
};

vecInv sphere(vec3 center, float radius) {
    return vecInv(
        center / radius,
        vec2(
            0.5 / radius,
            0.5 * (dot(center, center) / radius - radius)
        )
    );
}

// --- uniforms ---

// construction
const int SPHERE_MAX = 256;
uniform vecInv sphere_list[SPHERE_MAX];

// view
uniform vec2 resolution;
uniform float shortdim;

// controls
uniform vec2 ctrl;
uniform vec2 radius;
uniform float opacity;
uniform float highlight;
uniform int layer_threshold;
uniform bool use_test_construction;

// light and camera
const float focal_slope = 0.3;
const vec3 light_dir = normalize(vec3(2., 2., 1.));
const float ixn_threshold = 0.005;

// --- sRGB ---

// map colors from RGB space to sRGB space, as specified in the sRGB standard
// (IEC 61966-2-1:1999)
//
//   https://www.color.org/sRGB.pdf
//   https://www.color.org/chardata/rgb/srgb.xalter
//
// in RGB space, color value is proportional to light intensity, so linear
// color-vector interpolation corresponds to physical light mixing. in sRGB
// space, the color encoding used by many monitors, we use more of the value
// interval to represent low intensities, and less of the interval to represent
// high intensities. this improves color quantization

float sRGB(float t) {
    if (t <= 0.0031308) {
        return 12.92*t;
    } else {
        return 1.055*pow(t, 5./12.) - 0.055;
    }
}

vec3 sRGB(vec3 color) {
    return vec3(sRGB(color.r), sRGB(color.g), sRGB(color.b));
}

// --- shading ---

struct taggedFrag {
    int id;
    vec4 color;
    vec3 pt;
    vec3 normal;
};

taggedFrag[2] sort(taggedFrag a, taggedFrag b) {
    taggedFrag[2] result;
    if (a.pt.z > b.pt.z) {
        result[0] = a;
        result[1] = b;
    } else {
        result[0] = b;
        result[1] = a;
    }
    return result;
}

taggedFrag sphere_shading(vecInv v, vec3 pt, vec3 base_color, int id) {
    // the expression for normal needs to be checked. it's supposed to give the
    // negative gradient of the lorentz product between the impact point vector
    // and the sphere vector with respect to the coordinates of the impact
    // point. i calculated it in my head and decided that the result looked good
    // enough for now
    vec3 normal = normalize(-v.sp + 2.*v.lt.s*pt);
    
    float incidence = dot(normal, light_dir);
    float illum = mix(0.4, 1.0, max(incidence, 0.0));
    return taggedFrag(id, vec4(illum * base_color, opacity), pt, normal);
}

// --- ray-casting ---

vec2 sphere_cast(vecInv v, vec3 dir) {
    float a = -v.lt.s * dot(dir, dir);
    float b = dot(v.sp, dir);
    float c = -v.lt.t;
    
    float scale = -b/(2.*a);
    float adjust = 4.*a*c/(b*b);
    
    if (adjust < 1.) {
        float offset = sqrt(1. - adjust);
        return vec2(
            scale * (1. - offset),
            scale * (1. + offset)
        );
    } else {
        // these parameters describe points behind the camera, so the
        // corresponding fragments won't be drawn
        return vec2(-1., -1.);
    }
}

void main() {
    const int sphere_cnt = 3;
    
    vec2 scr = (2.*gl_FragCoord.xy - resolution) / shortdim;
    vec3 dir = vec3(focal_slope * scr, -1.);
    
    // initialize two spheres
    vecInv sphere_list_internal [sphere_cnt];
    if (use_test_construction) {
        /* DEBUG */
        // spheres 0 and 1 are identical in the test construction hard-coded
        // here and the construction passed in through uniforms. sphere 2 has
        // a different radius in the construction; we can use that to show that
        // the switch is working
        sphere_list_internal[0] = sphere(vec3(0.5, 0.5, -5. + ctrl.x), radius.x);
        sphere_list_internal[1] = sphere(vec3(-0.5, -0.5, -5. + ctrl.y), radius.y);
        sphere_list_internal[2] = sphere(vec3(-0.5, 0.5, -5.), 0.5);
    } else {
        for (int i = 0; i < sphere_cnt; ++i) {
            sphere_list_internal[i] = sphere_list[i];
        }
    }
    vec3 color_list [sphere_cnt];
    color_list[0] = vec3(1., 0.25, 0.);
    color_list[1] = vec3(0., 0.25, 1.);
    color_list[2] = vec3(0.25, 0., 1.0);
    
    // cast rays through the spheres
    vec2 depth_pairs [sphere_cnt];
    taggedFrag frags [2*sphere_cnt];
    int frag_cnt = 0;
    for (int i = 0; i < sphere_cnt; ++i) {
        vec2 hit_depths = sphere_cast(sphere_list_internal[i], dir);
        if (!isnan(hit_depths[0])) {
            frags[frag_cnt] = sphere_shading(sphere_list_internal[i], hit_depths[0] * dir, color_list[i], i);
            ++frag_cnt;
        }
        if (!isnan(hit_depths[1])) {
            frags[frag_cnt] = sphere_shading(sphere_list_internal[i], hit_depths[1] * dir, color_list[i], i);
            ++frag_cnt;
        }
    }
    
    // sort the fragments by depth, using an insertion sort
    for (int take = 1; take < frag_cnt; ++take) {
        taggedFrag pulled = frags[take];
        for (int put = take; put >= 0; --put) {
            if (put < 1 || frags[put-1].pt.z >= pulled.pt.z) {
                frags[put] = pulled;
                break;
            } else {
                frags[put] = frags[put-1];
            }
        }
    }
    
    // highlight intersections and cusps
    for (int i = frag_cnt-1; i >= 1; --i) {
        // intersections
        taggedFrag frag0 = frags[i];
        taggedFrag frag1 = frags[i-1];
        float ixn_sin = length(cross(frag0.normal, frag1.normal));
        vec3 disp = frag0.pt - frag1.pt;
        float ixn_dist = max(
            abs(dot(frag1.normal, disp)),
            abs(dot(frag0.normal, disp))
        ) / ixn_sin;
        float ixn_highlight = 0.5 * highlight * (1. - smoothstep(2./3.*ixn_threshold, 1.5*ixn_threshold, ixn_dist));
        frags[i].color = mix(frags[i].color, vec4(1.), ixn_highlight);
        frags[i-1].color = mix(frags[i-1].color, vec4(1.), ixn_highlight);
        
        // cusps
        float cusp_cos = abs(dot(dir, frag0.normal));
        float cusp_threshold = 2.*sqrt(ixn_threshold * sphere_list_internal[frag0.id].lt.s);
        float cusp_highlight = highlight * (1. - smoothstep(2./3.*cusp_threshold, 1.5*cusp_threshold, cusp_cos));
        frags[i].color = mix(frags[i].color, vec4(1.), cusp_highlight);
    }
    
    // composite the sphere fragments
    vec3 color = vec3(0.);
    for (int i = frag_cnt-1; i >= layer_threshold; --i) {
        if (frags[i].pt.z < 0.) {
            vec4 frag_color = frags[i].color;
            color = mix(color, frag_color.rgb, frag_color.a);
        }
    }
    outColor = vec4(sRGB(color), 1.);
}