// https://www.shadertoy.com/view/NsjSDK

/// Settings

// Use a fixed-point calculations while ray marching
// to prevent cumulative floating-point addition error.
// I originally implemented this to fix an artifact on the shadow side of the cube,
// and it worked, but I appear to have accidentally fixed the issue
// with the floating-point version too since then.
// Neither one appears to be more performant at this point,
// but I'd lean towards using the fixed point version to avoid any more
// of the aforementioned artifacts.
#define FIXED_POINT_MARCH
// 4x supersampling
#define SUPERSAMPLING
// All perspective introduces some form of distortion of the image.
// Non-curvilinear perspective preserves x/y parallel lines,
// whereas curvilinear perspective preserves the distances to objects.
// I prefer curvilinear perspective when the camera is zoomed in
// to minimize distortion.
// If you set the FOV to very high (up to 180 degrees), the scene will be
// cut off in a circle, just as a byproduct of how curvilinear perspective works.
// (It projects the surface of half the surface of a sphere around the camera onto a circle,
// whereas linear perspective projects onto a plane in front of the camera.)
#define CURVILINEAR
// Very large view distance, making the background of the scene much nicer,
// and preventing the shadows from getting cut off.
// However, rays sailing off into infinity is very slow
// (because the ground plane is infinite, the ray march distance interval
// doesn't increase very quickly or even shrinks, meaning each ray
// that's not substantially above the horizon takes nearly MAX_STEPS),
// so this option is not recommended in combination with supersampling
// unless you're rendering a still image.
//#define LONG_VIEW
// The width of the field of view. (180 degrees / FOV)
// Larger numbers suffer from less distortion at the edges of the image,
// but make the scene appear more zoomed-in.
#define FOV 3.2
// The distance of the camera from the scene.
// The ball sails around in a circle with a distance of 5 and radius 4,
// so increasing this distance will result in the ball passing through the camera.
#define CAMERA_DISTANCE 3.99


// The minimum distance resolution used by the fixed-point ray marching algorithm
// is 1/ this number. You can get away with as low as 1024 without much noticable difference.
// You'd think a very high resolution would mean that the algorithm runs slower,
// but in practice the distance to the object halves each iteration, so
// the number of iterations is more-or-less logarithmic with regards to this number.
// The only place it's slower is when a ray passes next to the edges of an object,
// but that's what the resolution is *for*: preventing graphical artifacts at the edges of objects.
#define FIXED_POINT_RESOLUTION 4096
#define FIXED_POINT_RESOLUTION_F 4096.
// The minimum distance from a surface that counts as a collision by the ray marching algorithm.
#ifdef FIXED_POINT_MARCH
#define MIN_DISTANCE (1./FIXED_POINT_RESOLUTION_F)
#else
// Any value much lower than this breaks the algorithm due to rounding error.
#define MIN_DISTANCE 0.001
#endif

#ifdef LONG_VIEW
// The greatest view distance that I can tell makes any difference at 800x450.
// My GPU runs this at 60 FPS with supersampling disabled at 800x450.
#define MAX_STEPS 3000
// ^ the maximum number of iterations taken by ray marching
#define MAX_DISTANCE 400.
#define FIXED_POINT_MAX_DISTANCE (400 * FIXED_POINT_RESOLUTION)
#else
// The shortest view distance which I don't think cuts off *too* much of the stage.
// The ideal would be to get my GPU to run this at 60 FPS at 1920x1080 with supersampling,
// but so far I haven't been able to optimize it to that point,
// and I'm not sure I even *could*, with my GPU and as slow of an algorithm as ray marching.
#define MAX_STEPS 200
#define MAX_DISTANCE 20.
#define FIXED_POINT_MAX_DISTANCE (20 * FIXED_POINT_RESOLUTION)
#endif

#define PI 3.14159

// problems to figure out: raymarching performance, numerical precision, math

// sRGB to linear RGB and vice versa
vec3 from_srgb(vec3 srgb) {
    // approximation
    return vec3(pow(srgb.x, 1./2.2), pow(srgb.y, 1./2.2), pow(srgb.z, 1./2.2));
}

vec3 from_linear(vec3 linear_rgb)
{
    // copied from somewhere on the internet
    bvec3 cutoff = lessThan(linear_rgb, vec3(0.0031308));
    vec3 higher = vec3(1.055)*pow(linear_rgb, vec3(1.0/2.4)) - vec3(0.055);
    vec3 lower = linear_rgb * vec3(12.92);

    return mix(higher, lower, cutoff);
}

vec3 to_polar(vec3 cart) {
    float r = length(cart);
    return vec3(r, acos(cart.y / r), atan(cart.z, cart.x));
}

vec3 from_polar(vec3 pol) {
    float xz = sin(pol.y);
    return pol.x * vec3(cos(pol.z) * xz, cos(pol.y), sin(pol.z) * xz);
}

vec3 realign(vec3 axis, vec3 vec) {
    vec3 axis_p = vec3(0., to_polar(axis).yz);
    vec3 vec_p = to_polar(vec);
    return from_polar(vec_p - axis_p);
}

// Signed distance functions of various shapes.
// These are fundamental to how ray marching works,
// but I will not be explaining how they (or ray marching) works any further.
float sdf_segment(vec3 s1, vec3 s2, vec3 p) {
    // Translate relative to s1.
    vec3 s1s2 = s2 - s1;
    vec3 s1p = p - s1;
    vec3 s_direction = normalize(s1s2);
    // Distance along s nearest to p.
    float h = dot(s_direction, s1p);
    // Point along s nearest to p.
    vec3 s1x = s_direction * clamp(h, 0., length(s1s2));
    // Distance from p to x.
    return distance(s1p, s1x);
}

float sdf_plane(mat4 trans, vec3 pos) {
    vec4 p1w = trans * vec4(0., 0., 0., 1.);
    vec4 p2w = trans * vec4(1., 0., 0., 1.);
    vec4 p3w = trans * vec4(0., 0., 1., 1.);
    vec3 p1 = p1w.xyz / p1w.w;
    vec3 p2 = p2w.xyz / p2w.w;
    vec3 p3 = p3w.xyz / p3w.w;
    vec3 n = normalize(cross(p2 - p1, p3 - p1));
    vec3 rp = pos - p1;
    return dot(rp, n);
}

float sdf_pill(vec3 s1, vec3 s2, float radius, vec3 p) {
    return sdf_segment(s1, s2, p) - radius;
}

float sdf_sphere(vec3 sphere_pos, float sphere_radius, vec3 pos) {
    return distance(pos, sphere_pos) - sphere_radius;
}

float sdf_cube(
    mat4 trans,
    vec3 pos
) {
    // positive faces of cube
    mat4 f1 = trans * mat4(
        0.,  1., 0., 0.,
        1.,  0., 0., 0.,
        0.,  0., 1., 0.,
        0.5, 0., 0., 1.
    );
    // WTF: why does this y have to be negative?
    mat4 f2 = trans * mat4(
        1., 0., 0., 0.,
        0., 1., 0., 0.,
        0., 0., 1., 0.,
        0., -0.5, 0., 1.
    );
    mat4 f3 = trans * mat4(
        1., 0., 0., 0.,
        0., 0., 1., 0.,
        0., 1., 0., 0.,
        0., 0., 0.5, 1.
    );
    
    // signed distances from positive cube faces
    float df1 = sdf_plane(f1, pos);
    float df2 = sdf_plane(f2, pos);
    float df3 = sdf_plane(f3, pos);
    
    // distances between parallel cube faces
    float lpf1 = length(trans * vec4(1., 0., 0., 0.));
    float lpf2 = length(trans * vec4(0., 1., 0., 0.));
    float lpf3 = length(trans * vec4(0., 0., 1., 0.));
    
    // signed distances from both parallel cube faces
    // * if positive, positive distance from plane
    // * if negative:
    //   * remeasure difference from opposite face (which is +1 away)
    //   * if distance from opposite face is positive, then we are in the cube
    //   * otherwise we are outside the cube
    
    float dpf1 = max(df1, -(df1 + lpf1));
    float dpf2 = max(df2, -(df2 + lpf2));
    float dpf3 = max(df3, -(df3 + lpf3));
    
    vec3 dpfs = vec3(dpf1, dpf2, dpf3);
    
    // if the point is outside the cube, then a point within two faces adds 0 distance;
    // if the point is inside the cube, then a point within the faces does add distance.
    if (all(lessThan(dpfs, vec3(0.)))) {
        // the point is inside the cube.
        return -length(dpfs);
    }
    
    return length(max(dpfs, 0.));
}

float df_cube_minus_sphere(mat4 trans, float sphere_radius, vec3 pos) {
    // TODO: SDF for this object
    float sdcube = sdf_cube(trans, pos);
    if (sdcube >= MIN_DISTANCE) {
        return sdcube;
    }
    
    // we're inside the cube
    vec3 sphere_pos = (trans * vec4(0., 0., 0., 1.)).xyz;
    float sdsphere = sdf_sphere(sphere_pos, sphere_radius, pos);
    return max(sdcube, -sdsphere);
}

float distance_to_objects(out int material, vec3 pos) {
    pos.z -= CAMERA_DISTANCE;
    //pos = realign(vec3(1., 0., 0.), pos);
    //pos = from_polar(to_polar(pos));
    
    vec3 sphere_pos = 5.0 * vec3(sin(iTime), 0.0, cos(iTime));
    float sphere_radius = 1.0;
    float sphere_dist = abs(sdf_sphere(sphere_pos, 1.0, pos));
    
    mat4 ground_trans = mat4(
        1.,  0., 0., 0.,
        0.,  1., 0., 0.,
        0.,  0., 1., 0.,
        0., -1., 0., 1.
    );
    float ground_dist = abs(sdf_plane(ground_trans, pos));
    
    vec3 pill_pos1 = vec3(-1.0, -0.5, 2.7);
    vec3 pill_pos2 = vec3(0.4, -0.9, 2.5);
    float pill_radius = 0.1;
    float pill_dist = abs(sdf_pill(pill_pos1, pill_pos2, pill_radius, pos));
        
    mat4 cube_trans = mat4(
        0.6, 0., 0., 0.,
        0., 0.6, 0., 0.,
        0., 0., 0.6, 0.,
        0.2 - 1.3 * sin(iTime/2.), 0.3 - 0.3 * cos(iTime/3.), 2.5 + 0.2 * sin((iTime+0.5)*2.), 1.
    );
    float cube_angle = iTime;
    mat4 cube_rotate = mat4(
        cos(cube_angle), -sin(cube_angle), 0., 0.,
        sin(cube_angle), cos(cube_angle), 0., 0.,
        0., 0., 1., 0.,
        0., 0., 0., 1.
    );
    //float cube_dist = abs(sdf_cube(cube_trans, pos));
    float cube_dist = df_cube_minus_sphere(cube_trans * cube_rotate, 0.39, pos);
    
    // HACK: What function generalizes this pattern?
    // Obviously min for the distances, but what about the materials?
    if (sphere_dist < ground_dist && sphere_dist < pill_dist && sphere_dist < cube_dist) {
        material = 0;
        return sphere_dist;
    }
    if (ground_dist < pill_dist && ground_dist < cube_dist) {
        material = 1;
        return ground_dist;
    }
    if (pill_dist < cube_dist) {
        material = 2;
        return pill_dist;
    }
    material = 3;
    return cube_dist;
}

vec3 material_color(int mat) {
    if (mat == 0) {
        return from_srgb(vec3(0.0, 0.19, 0.56));
    }
    if (mat == 1) {
        return vec3(1.0, 1.0, 1.0);
    }
    if (mat == 2) {
        return from_srgb(vec3(0.50, 1.00, 0.40));
    }
    if (mat == 3) {
        return from_srgb(vec3(0.83, 0.69, 0.22));
    }
    return vec3(1.0, 0.0, 1.0);
}

#ifndef FIXED_POINT_MARCH
float ray_march(out int mat, vec3 ray_origin, vec3 ray_direction) {
    vec3 position_along_ray = ray_origin;
    float dist = 1./0.;
    mat = -1;
    for (int i = 0; i < MAX_STEPS && distance(ray_origin, position_along_ray) < MAX_DISTANCE && dist > MIN_DISTANCE; i++) {
        dist = distance_to_objects(mat, position_along_ray);
        position_along_ray += ray_direction * dist;
    }
    if (distance(ray_origin, position_along_ray) > MAX_DISTANCE) {
        return 1.0/0.0;
    }
    return distance(ray_origin, position_along_ray);
}
#else
float ray_march(out int mat, vec3 ray_origin, vec3 ray_direction) {
    int travel_dist = 0;
    int delta = 1;
    for (int i = 0; i < MAX_STEPS && delta > 0; i++) {
        float dist = distance_to_objects(mat, ray_origin + ray_direction * float(travel_dist) / FIXED_POINT_RESOLUTION_F);
        delta = min(int(floor(dist * FIXED_POINT_RESOLUTION_F)), FIXED_POINT_MAX_DISTANCE - travel_dist);
        travel_dist += delta;
    }
    if (travel_dist >= FIXED_POINT_MAX_DISTANCE) {
        return 1./0.;
    }
    return float(travel_dist) / FIXED_POINT_RESOLUTION_F;
}
#endif

vec3 normal(vec3 pos) {
    // Magic number determined by tinkering with it until stuff worked.
    vec2 delta = vec2(0.0001, 0.);
    int _mat = -1;
    float dist = distance_to_objects(_mat, pos);
    vec3 dq = (dist - vec3(
        distance_to_objects(_mat, pos - delta.xyy),
        distance_to_objects(_mat, pos - delta.yxy),
        distance_to_objects(_mat, pos - delta.yyx)
    )) / delta.x;
    
    return normalize(dq);
}

float light(vec3 pos) {
    vec3 source = vec3(0.5, 1.5, -1.0);
    float intensity = dot(normalize(source - pos), normal(pos));
    vec3 direction = normalize(source - pos);
    int _mat = -1;
    // 10 is a magic number determined by experimentation.
    // Anything less generates noticable artifacts due to rounding error (depending on on MIN_DISTANCE).
    if (ray_march(_mat, pos + direction * (10. * MIN_DISTANCE), direction) < distance(source, pos)) {
        intensity *= 0.1;
    }
    if (intensity < 0.005 && length(pos) < 1./0.) {
        return 0.005;
    }
    return clamp(intensity, 0., 1.);
}

vec3 project_spherical(vec2 pos) {
    return vec3(pos, sqrt(1.0 - pow(pos.x, 2.0) - pow(pos.y, 2.0)));
}

vec3 project_planar(vec2 pos) {
    return normalize(vec3(pos.x, pos.y, 1.0));
}

// Cast a ray in a direction and figure out the color of the material.
vec3 color_ray(vec3 dir) {
    int mat = -1;
    vec3 pos = vec3(0.);
    vec3 color = vec3(0.);
    //for (int i = 1; i < 3; i++) {
    int i = 1;
        float dist = ray_march(mat, vec3(0.0), dir);
        pos += dir * dist;
        color += material_color(mat) * light(pos) / float(i);
    //    dir = reflect(dir, normal(pos));
    //    pos += 10. * MIN_DISTANCE * dir;
    //}
    return color;
}

vec3 color_frame(vec2 loc) {
    #ifdef CURVILINEAR
    vec3 dir = project_spherical(loc);
    #else
    vec3 dir = project_planar(loc);
    #endif
    return color_ray(dir);
}

// 4x the samples means 4x better-looking edges and 4x slower.
vec3 supersample(vec2 loc, float increment) {
    increment /= 2. * FOV;
    vec2 d = vec2(increment, -increment);
    // As far as I can tell, antialiasing algorithms usually put the samples at strange locations,
    // but I neither know why nor what the actual offsets are.
    return (color_frame(loc + d.xx) + color_frame(loc + d.xy) + color_frame(loc + d.yx) + color_frame(loc + d.yy)) / 4.;
}

void mainImage(out vec4 frag_color, in vec2 frag_coord) {
    // Scale from pixel coordinates to scene coordinates.
    vec2 uv = (frag_coord / iResolution.xy * 2.0) - 1.0;
    vec2 loc;
    // Cut off the part of the scene which doesn't fit due to the aspect ratio.
    if (iResolution.x > iResolution.y) {
        loc = vec2(uv.x, uv.y * iResolution.y / iResolution.x);
    } else {
        loc = vec2(uv.x * iResolution.x / iResolution.y, uv.y);
    }
    
    vec2 xy = loc / FOV;

    #ifdef SUPERSAMPLING
    // FIXME: 1/iResolution.x is incorrect if the scene is larger than it is tall.
    vec3 color = supersample(xy, 1. / iResolution.x);
    #else
    vec3 color = color_frame(xy);
    #endif
    int mat = -1;
    //color = vec3(distance_to_objects(mat, 5. * project_spherical(xy)) / MAX_DISTANCE);
    //color = vec3(ray_march(mat, vec3(0.), project_spherical(xy)) / MAX_DISTANCE);
    vec3 dir = project_spherical(xy);
    float dist = ray_march(mat, vec3(0.), dir);
    //color = material_color(mat) * dist / MAX_DISTANCE;
    //color = (normal(dist * dir) + vec3(1.0)) / 2.;
    
    frag_color = vec4(from_linear(color), 0.0);
}