// References. Big thanks to the authors of these awesome shaders/tutorials.
// This has been a good exercise helping me reviewing basic concepts of 
// geometry transformation, raymarching and shading model.
// 
// RayMarching starting point by BigWIngs. https://shadertoy.com/view/WtGXDD
// Bending Light https://www.youtube.com/watch?v=NCpaaLkmXI8
// Sphere in a box https://www.shadertoy.com/view/3sycDm
// Ray Intersector https://iquilezles.org/articles/intersectors
// Creation by Silexars  https://shadertoy.com/view/XsXXDn
//
// Will not include in this excercise but will consider practicing later
// 1. Photo editing / image processing algorithms, grain
// 2. DOF and bokeh
// 3. Black Body Radiation -> Dispersion Function (Sellmeier approximation) -> Wavelength -> RGB

// Overall Graphics
#define AA 2
#define DIM_LEVEL 0.4
#define OBJ_AMP 1.0 

//RayIntersect
#define SURF_DIST .001
// The offset to work with interection on both side
#define OFFSET 15.
#define FLOAT_MAX 3.402823466e+38

// Bouncing and spectral setup
#define IOR 2.45
#define ABB 0.06
#define MATERIAL_DECAY 0.1
#define MAX_BOUNCES 10
#define SPECTRAL_BANDS 3

// Helper
#define S smoothstep
#define T iTime
#define PI 3.14159265358
#define keyToggle(ascii)  ( texelFetch(iChannel1,ivec2(ascii,2),0).x > 0.)


const vec3 SP_COL[SPECTRAL_BANDS] = vec3[](vec3(1.0, 0.0, 0.0), 
                                           vec3(0.0, 1.0, 0.0),
                                           vec3(0.0, 0.0, 1.0));

// stores intersection result
struct ObjInter
{
    bool inter;
    bool inside;
    float d;
    vec3 n;
};

float random(in vec3 p)
{
    return fract(sin(dot(p,vec3(.7434,.4234, .793)))*33348.541133);
}

float noise(in vec3 p)
{
	p.z += iTime * .4;
	
    vec3 i = floor(p);
	vec3 f = fract(p); 
	f *= f * (3.-2.*f);

    vec2 c = vec2(0,1);

    return mix(
		mix(mix(random(i + c.xxx), random(i + c.yxx),f.x),
			mix(random(i + c.xyx), random(i + c.yyx),f.x),
			f.y),
		mix(mix(random(i + c.xxy), random(i + c.yxy),f.x),
			mix(random(i + c.xyy), random(i + c.yyy),f.x),
			f.y),
		f.z);
}

mat2 Rot(float a) {
    float s=sin(a), c=cos(a);
    return mat2(c, -s, s, c);
}

mat3 Rot3D(vec3 ang)
{
    mat3 mx = mat3(
			1.0,		0.0,		0.0,
			0.0,		cos(ang.x),	-sin(ang.x),
			0.0,		sin(ang.x),	cos(ang.x));
    mat3 my = mat3(
			cos(ang.y), 0.0,		sin(ang.y),
			0.0,		1.0,		0.0,
			-sin(ang.y),0.0,		cos(ang.y));
    mat3 mz = mat3(
			cos(ang.z), -sin(ang.z),0.0,
			sin(ang.z),	cos(ang.z),	0.0,
			0.0,		0.0,		1.0);
        
    return mz*my*mx;
}

// Intersector https://iquilezles.org/articles/intersectors
float sphIntersect( in vec3 ro, in vec3 rd, in vec4 sph )
{
	vec3 oc = ro - sph.xyz;
	float b = dot( oc, rd );
	float c = dot( oc, oc ) - sph.w*sph.w;
	float h = b*b - c;
	if( h<0.0 ) return -1.0;
	return -b - sqrt( h );
}

vec3 sphNormal( in vec3 pos, in vec4 sph )
{
    return normalize(pos-sph.xyz);
}

float roundedboxIntersect( in vec3 ro, in vec3 rd, in vec3 size, in float rad )
{
    // bounding box
    vec3 m = 1.0/rd;
    vec3 n = m * ro;
    vec3 k = abs(m)*(size+rad);
    vec3 t1 = -n - k;
    vec3 t2 = -n + k;
    float tN = max( max( t1.x, t1.y ), t1.z );
    float tF = min( min( t2.x, t2.y ), t2.z );
    if( tN>tF || tF<0.0) return -1.0;
    float t = tN;

    // convert to first octant
    vec3 pos = ro+t*rd;
    vec3 s = sign(pos);
    ro  *= s;
    rd  *= s;
    pos *= s;
        
    // faces
    pos -= size;
    pos = max( pos.xyz, pos.yzx );
    if( min(min(pos.x,pos.y),pos.z) < 0.0 ) return t;

    // some precomputation
    vec3 oc = ro - size;
    vec3 dd = rd*rd;
    vec3 oo = oc*oc;
    vec3 od = oc*rd;
    float ra2 = rad*rad;

    t = 1e30;        

    // corner
    {
    float b = od.x + od.y + od.z;
    float c = oo.x + oo.y + oo.z - ra2;
    float h = b*b - c;
    if( h>0.0 ) t = -b-sqrt(h);
    }
    // edge X
    {
    float a = dd.y + dd.z;
    float b = od.y + od.z;
    float c = oo.y + oo.z - ra2;
    float h = b*b - a*c;
    if( h>0.0 )
    {
        h = (-b-sqrt(h))/a;
        if( h>0.0 && h<t && abs(ro.x+rd.x*h)<size.x ) t = h;
    }
    }
    // edge Y
    {
    float a = dd.z + dd.x;
    float b = od.z + od.x;
    float c = oo.z + oo.x - ra2;
    float h = b*b - a*c;
    if( h>0.0 )
    {
        h = (-b-sqrt(h))/a;
        if( h>0.0 && h<t && abs(ro.y+rd.y*h)<size.y ) t = h;
    }
    }
    // edge Z
    {
    float a = dd.x + dd.y;
    float b = od.x + od.y;
    float c = oo.x + oo.y - ra2;
    float h = b*b - a*c;
    if( h>0.0 )
    {
        h = (-b-sqrt(h))/a;
        if( h>0.0 && h<t && abs(ro.z+rd.z*h)<size.z ) t = h;
    }
    }

    if( t>1e19 ) t=-1.0;
    
    return t;
}

// normal of a rounded box
vec3 roundedboxNormal( in vec3 pos, in vec3 siz, in float rad )
{
    return sign(pos)*normalize(max(abs(pos)-siz,0.0));
}

ObjInter roundedBoxIntersect(in vec3 roB, in vec3 roF, in vec3 rd, vec3 box_c, vec3 size, float r, mat3 m) {
    ObjInter box;
    box.inside = false;
    
    roB -= box_c;
    roF -= box_c;
    
    roB = m * roB;
    roF = m * roF;
    rd = m * rd;
    
    vec3 p;
    float d1 = roundedboxIntersect(roB, rd, size, r);
    float d2 = roundedboxIntersect(roF, -rd, size, r);
    
    if (d1 > OFFSET) {
        // outside
        p = roB + rd*d1;
        box.d = d1 - OFFSET;
        box.n = roundedboxNormal(p, size, r) * m;
        box.inter = true;
    } else if (d2 > 0.0 && d2 < OFFSET) {
        // inside
        p = roF - rd*d2;
        box.d = OFFSET - d2;
        box.n = -roundedboxNormal(p, size, r) * m;
        box.inside = true;
        box.inter = true;
    } else {
        box.inter = false;
    }
    return box;
}

ObjInter sphereIntersect(in vec3 roB, in vec3 roF, in vec3 rd, vec4 para) {
    ObjInter sph;
    sph.inside = false;
    
    vec3 p;
    float d1 = sphIntersect(roB, rd, para);
    float d2 = sphIntersect(roF, -rd, para);
    
    if (d1 > OFFSET) {
        // outside
        p = roB + rd*d1;
        sph.d = d1 - OFFSET;
        sph.n = sphNormal(p, para);
        sph.inter = true;
        
    } else if (d2 > 0.0 && d2 < OFFSET) {
        // d1 < d2
        // inside
        p = roF - rd*d2;
        sph.d = OFFSET - d2;
        sph.n = -sphNormal(p, para);
        sph.inside = true;
        sph.inter = true;
    } else {
        sph.inter = false;
    }
    return sph;
}

float planeIntersect(vec3 ro, vec3 rd, vec3 nor, float dist) {
  float denom = dot(rd, nor);
  float t = -(dot(ro, nor) + dist) / denom;

  return t;
}

// Source: https://iquilezles.org/articles/palettes
vec3 palette( in float t, in vec3 a, in vec3 b, in vec3 c, in vec3 d ) {
    return a + b*cos( 2.*PI*(c*t+d) );
}

vec3 paletteBackground(vec3 ro, vec3 rd) {
    vec3 col = palette(
        clamp(rd.z * 0.5 + 0.5, -1.0, 1.0),
        vec3(0.5, 0.5, 0.5),
        vec3(0.5, 0.5, 0.5),
        vec3(1.0, 1.0, 1.0),
        vec3(0.0, 0.1, 0.20));
  return col;
}

// Creation by Silexars  https://shadertoy.com/view/XsXXDn
vec3 litPlane(vec3 ro, vec3 rd) {
    vec3 col;
    
    ro.xz *= Rot(-0.1*T);
    rd.xz *= Rot(-0.1*T);

    // get the plane the BG is displaying on
    float d = planeIntersect(ro, rd, normalize(vec3(0, -1, 0)), 20.);

    if (d > 0.0) {
        vec3 pp = ro + rd * d;
        
        /* position */
        vec2 f = pp.xz;
        vec2 U;
        vec2 r = vec2(40, 40);
        vec2 p = f/r;

        p.x*= r.x/r.y;

        float t = 0.2*T, l = length(p);
        
        
        for( int i=0; i<3; i++ ) {
            t +=.02;
            U = f/r + normalize(p) * ( sin(t)+1.) * abs(sin(l*9.-t*2.));
            col[i] = .01 / length( abs(fract(U)-.5) );
        }
        
        col = col/l;
    }
    return col * smoothstep(40., 30., d);
}

vec3 getBackgroundCol(in vec3 ro, in vec3 rd) {
    vec3 col;

    // ig's Palette
    col = paletteBackground(ro, rd);
    
    // plane
    col += litPlane(ro, rd) * 0.5;
    
    // revert gamma correction since we'll do final correction
    return pow(col, vec3(2.2));
}

// RayMarching is too slow
// Use intersector instead
// now the inside para is used to determine where is the start point located
bool RayMarchIntersect(in vec3 ro, in vec3 rd, in float side, out float d, out vec3 oNor) {
    // Set up the scene
    const int max_obj = 11;
    ObjInter objs[max_obj];
    
    // Offset on both direction
    // Note: only work accurately on simple geometry
    // For shape like a torus, even sometime a ray has
    // four intersection, this only give you two of them
    
    // move backwards, look forwards
    vec3 roB = ro - rd*OFFSET;
    // move forwards, look backwards
    vec3 roF = ro + rd*OFFSET;
    
    // Rounded Box 0
    float box_r = 0.02;
    vec3 box_c = vec3(0., 0., 0.);
    vec3 box_size = vec3(4, 4, 4);
    mat3 box_m = Rot3D(vec3(0, 0., 0.));
    ObjInter box = roundedBoxIntersect(roB, roF, rd, box_c, box_size, box_r, box_m);
    //objs[0] = box;
    

    // Rounded Box 1
    float box_r1 = .02;
    vec3 box_c1 = vec3(0, 0, 0);
    vec3 box_size1 = vec3(0.9);
    mat3 box_m1 = Rot3D(vec3(0, 0.1*T, 0.));
    //mat3 box_m = Rot3D(vec3(atan(1./sqrt(2.)), 0, 1./4.*PI));
    ObjInter box1 = roundedBoxIntersect(roB, roF, rd, box_c1, box_size1, box_r1, box_m1);
    objs[1] = box1;
    
    // Sphere 1
    vec3 loc = vec3(0., 1.1, 0.);
    loc.xy = Rot(0.2*T) * loc.xy;
    vec4 sph_para = vec4(loc, 1.1);
    ObjInter sph = sphereIntersect(roB, roF, rd, sph_para);
    objs[2] = sph;
    
    // Sphere 2
    vec3 loc2 = vec3(0., -1.1, 0.);
    loc2.xy = Rot(0.2*T) * loc2.xy;
    vec4 sph_para2 = vec4(loc2, 1.1);
    ObjInter sph2 = sphereIntersect(roB, roF, rd, sph_para2);;
    objs[3] = sph2;

    // Sphere shell
    vec4 sph_para3 = vec4(0.0, 0.0, 0.0, 3.0);
    ObjInter sph3 = sphereIntersect(roB, roF, rd, sph_para3);;
    objs[4] = sph3;
    
    // a generic intersection handling
    // find the nearest intersection
    int inter_count = 0;
    ObjInter interObjs[max_obj];
    
    for (int i = 0; i < max_obj; i++) {
        if (objs[i].inter) {
            interObjs[inter_count] = objs[i];
            inter_count++;
        }
    }
    
    if (inter_count < 1) {
        return false;
    }
    
    int min_idx;
    float min_d = FLOAT_MAX;
    for (int i = 0; i < inter_count; i++) {
        if (interObjs[i].d < min_d) {
            min_idx = i;
            min_d = interObjs[i].d;
        }
    }
    
    d = interObjs[min_idx].d;
    oNor = interObjs[min_idx].n;
    
    return true;
}


// p = point
// l = lookat
// z = zoom
// r_para is the vector implicitly determine rotation and
// handedness of coordinate system.
vec3 GetRayDir(vec2 uv, vec3 p, vec3 l, vec3 r_para, float z) {
    vec3 f = normalize(l-p),
        // this intermediate vec r is to make sure u is in the correct plane
        r = normalize(cross(r_para, f)),
        u = cross(f,r),
        c = f*z,
        i = c + uv.x*r + uv.y*u,
        d = normalize(i);
    return d;
}

// p: first contact point on surface
// rd: subsurface rd
vec3 renderSub(in vec3 p, in vec3 n, in vec3 rd, in float ior) {
    
    // the first refraction entering the solid
    vec3 rdStart = rd;
    // offset in a little bit to 
    // determin inner ray start point
    vec3 pStart = p - n*SURF_DIST;
    
    // start with inner side
    float side = -1.;
    
    // to avoid multiple declaration
    float dStart;
    vec3 rdNext, pNext, nNext;
    
    // ray path length
    float innerLength = 0.;
    
    for (int i = 0; i < MAX_BOUNCES; i++) {
        
        // subsurface bounce, scale down the step
        
        //dStart = RayMarch(pStart, rdStart, side, .5);
        float dStart;
        vec3 nNext;
        bool inter = RayMarchIntersect(pStart, rdStart, side, dStart, nNext);
        
        if (side < 0.) {
            innerLength += dStart;
        }   

        // if ray leaves solid
        if (!inter) {
            break;
        }

        pNext = pStart + rdStart * dStart;
        
        // for simplicity, we consider the ray either
        // relfect OR refract, not both.
        rdNext = refract(rdStart, nNext, ior);
        //
        // check if rdOut.xyz are all 0
        // length(rdOut) also does the trick,
        // but slower than dot(v, v)
        if (dot(rdNext, rdNext) == 0.) {
            pStart = pNext + nNext*SURF_DIST;
            rdStart = reflect(rdStart, nNext);
        } else {
            pStart = pNext - nNext*SURF_DIST;
            rdStart = rdNext;
            side = -side;
            ior = 1. / ior;
        }
    }
    
    // dIn is the distance/length of internal ray
    float optDist = exp(-innerLength * MATERIAL_DECAY);
    
    return getBackgroundCol(pStart, rdStart) * optDist;
}

vec3 vignette(vec3 color, vec2 q, float v)
{
    color *= 0.3 + 0.8 * pow(16.0 * q.x * q.y * (1.0 - q.x) * (1.0 - q.y), v);
    return color;
}

vec3 desaturate(in vec3 c, in float a)
{
    float l = dot(c, vec3(1. / 3.));
    return mix(c, vec3(l), a);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // the same offset and scale for uv when rendering,
    // to make the center (0,0), and space [-.5 .5];
    vec2 mouse = (iMouse.xy-.5*iResolution.xy)/iResolution.y;
    
    // set camera to enough far so it
    // looks like isometric
    // need more work on a real transformation
    vec3 ro = vec3(0, -10., 0);
    vec3 lookat = vec3(0, 0, 0);
    vec3 r_para = vec3(0, 0, 1.);
    float zoom = 1.; 
    
    // space key will lock scence to preset
    
    // Thanks davidar@ for the fix on blank screen when mouse is at (0,0)
    if (keyToggle(32) || iMouse.xy == vec2(0)) {
        // reset and lock mouse to center
        mouse = vec2(0.0);
    }
    
    //ro.xy *= Rot(m.x*PI*1.);
    //ro.yz *= Rot(m.y*PI*1.);
    //ro *= Rot3D
    ro = Rot3D(vec3(-mouse.y*PI, 0.0, -mouse.x*PI)) * ro;
    
    vec3 tot = vec3(.0);
    for (int m = 0; m < AA; m++) {
        for (int n = 0; n < AA; n++) {
            // AA sampling offset
            vec2 o = vec2(float(m),float(n)) / float(AA) - 0.5;
            vec2 uv = ((fragCoord+o)-.5*iResolution.xy)/iResolution.y;

            vec3 rd = GetRayDir(uv, ro, lookat, r_para, zoom);

            // the cube map is a 2d collection of rd endpoints
            vec3 col;

            // To save computation time, we only march once 
            // when hitting the object the first time
            //float d = RayMarch(ro, rd, 1., 1.);
            float d;
            vec3 oNor;
            bool inter = RayMarchIntersect(ro, rd, 1.0, d, oNor);


            if (inter) {
                vec3 p = ro + rd * d;
                //vec3 n = GetNormal(p);
                vec3 refl = reflect(rd, oNor);

                // reflection
                vec3 reflOutside = getBackgroundCol(p, refl);

                vec3 refrCol;

                // rgb chromatic abbreviation
                for (int i = 0; i < SPECTRAL_BANDS; i++) {
                    float ior = mix(IOR-ABB, IOR+ABB, float(i)/float(SPECTRAL_BANDS-1));
                    vec3 rdSub = refract(rd, oNor, 1./ior);
                    refrCol += renderSub(p, oNor, rdSub, ior) * SP_COL[i];
                }

                // flipped fresnel, could actually revert
                // doesn't matter.
                float fresnel = pow(1.-dot(-rd, oNor), 4.);

                col = mix(refrCol, reflOutside, fresnel);
                col *= OBJ_AMP;
                //col = col*dif;
            } else {
                // dim background
                col = getBackgroundCol(ro, rd);
            }
            tot += col;
        }
    }
    
    tot /= float(AA*AA);
    
    // gamma correction
    tot = pow(tot, vec3(.4545));
    //col = pow(col, vec3(1./2.2))
    
    //post processing
    tot = desaturate(tot, 0.1);
    tot = vignette(tot, fragCoord / iResolution.xy, 0.2);
    
    fragColor = vec4(tot,1.0);
}