The Smooth Cubic Thingie

#version 300 es
// ╔═════════════╦════════════════╗
// ║ Marco Gomez ║ https://mgz.me ║
// ╚═════════════╩════════════════╝
precision highp float;

uniform vec2 resolution;
uniform vec2 mouselerp;
uniform float time;

out vec4 fragColor;

#define marchSteps 128
#define supersample 0
#define maxDist 15.0
#define surfDist 0.0001

const float PI = acos(-1.0);
const float TAU = PI * 2.0;

float osc(float s, float e, float t, float ts) {
  return (e - s) / 2.0 + s + sin(t * ts) * (e - s) * 0.5;
}

mat3 xrot(float t) {
  return mat3(1.0, 0.0, 0.0, 0.0, cos(t), -sin(t), 0.0, sin(t), cos(t));
}

mat3 yrot(float t) {
  return mat3(cos(t), 0.0, -sin(t), 0.0, 1.0, 0.0, sin(t), 0.0, cos(t));
}

mat3 zrot(float t) {
  return mat3(cos(t), -sin(t), 0.0, sin(t), cos(t), 0.0, 0.0, 0.0, 1.0);
}

float sdBox(vec3 p, vec3 b) {
  vec3 d = abs(p) - b;
  return (min(max(d.x, max(d.y, d.z)), 0.0) + length(max(d, 0.0)));
}

float planeDistance(vec3 pos) {
  float o = clamp(osc(-2.5, 2.5, time, 0.12), -1.25, 0.0);
  vec3 origin = vec3(0.0, 0.5 + o, 0.0);
  vec3 normal = vec3(0.0, 1.0, 0.0);
  vec3 delta = pos - origin;
  float prod = dot(delta, normal);
  return prod;
}

float smoothUnion(float d1, float d2, float k) {
  float h = clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0);
  return mix(d2, d1, h) - k * h * (1.0 - h);
}

float smoothIntersection(float d1, float d2, float k) {
  float h = clamp(0.5 - 0.5 * (d2 - d1) / k, 0.0, 1.0);
  return mix(d2, d1, h) + k * h * (1.0 - h);
}

float mapDistance(vec3 pos) {
  float t = time * 0.75;
  vec3 rpos = (pos - vec3(0.0, 1.0, 0.0));
  rpos *= zrot(t) * xrot(PI * 0.25) * yrot(PI * 0.25 + t);
  float cube = sdBox(rpos, vec3(0.5));
  float acut = sdBox(rpos, vec3(1.00, 0.45, 0.45));
  float bcut = sdBox(rpos, vec3(0.45, 1.00, 0.45));
  float ccut = sdBox(rpos, vec3(0.45, 0.45, 1.00));
  float carve = smoothUnion(acut, smoothUnion(bcut, ccut, 0.025), 0.025);
  float x = smoothIntersection(-carve, cube, 0.0125);
  float dist = smoothUnion(x, planeDistance(pos), 1.25);
  float o = clamp(osc(-2.5, 2.5, time, 0.24), 0.0, 1.0);
  dist = mix(dist, min(dist, carve + 0.35), o);
  return dist;
}

vec3 normal(vec3 p) {
  vec3 o = vec3(0.01, 0.0, 0.0);
  vec3 n = vec3(0.0);
  n.x = mapDistance(p + o) - mapDistance(p - o);
  n.y = mapDistance(p + o.zxy) - mapDistance(p - o.zyx);
  n.z = mapDistance(p + o.yzx) - mapDistance(p - o.yzx);
  return normalize(n);
}

float rayMarch(vec3 ro, vec3 rd) {
  float dist = 0.0;
  for (int i = 0; i < marchSteps; ++i) {
    vec3 pos = ro + rd * dist;
    float d = mapDistance(pos);
    dist += d;
    if (dist > maxDist || d <= surfDist) { break; }
  }
  return dist;
}

vec2 lMarch(vec3 ro, vec3 rd) {
  float dist = 0.0;
  float md = 2.0;
  float lt = 0.0;
  for (int i = 0; i < 32; ++i) {
    vec3 pos = ro + rd * dist;
    float d = mapDistance(pos);
    md = min(md, 4.0 * d / dist);
    dist += min(d, 0.08);
  }
  return vec2(dist, md);
}

float light(vec3 world, vec3 sn, vec3 lpos) {
  vec3 ldel = world + sn * 0.01 - lpos;
  float ldist = length(ldel);
  ldel /= ldist;
  vec2 lt = lMarch(lpos, ldel);
  float lm = 1.0;
  if (lt.x < ldist) { lm = lt.y; }
  float lp = max(dot(ldel, -sn), 0.0);
  float fl = lp * lm / (1.0 + ldist * ldist * 0.1);
  return fl;
}

mat3 calcLookAtMatrix(vec3 origin, vec3 target, float roll) {
  vec3 rr = vec3(sin(roll), cos(roll), 0.0);
  vec3 ww = normalize(target - origin);
  vec3 uu = normalize(cross(ww, rr));
  vec3 vv = normalize(cross(uu, ww));
  return mat3(uu, vv, ww);
}

float getSoftShadow(vec3 ro, vec3 rd, float tmin, float tmax, float k) {
  float res = 1.0;
  float ph = 1e20;
  for(float t = tmin; t < tmax; ) {
    float h = mapDistance(ro + rd * t);
    if (h < 0.001) { return 0.0; }
    float y = h * h / (2.0 * ph);
    float d = sqrt(h * h - y * y);
    res = min(res, k * d / max(0.0, t-y));
    ph = h;
    t += h;
  }
  return res;
}

float getAmbientOcclusion(vec3 pos, vec3 nor) {
  float occ = 0.0;
  float sca = 1.0;
  for (int i = 0; i < 5; i++) {
    float h = 0.001 + 0.25 * float(i) / 4.0;
    float d = mapDistance(pos + h * nor);
    occ += (h - d) * sca;
    sca *= 0.98;
  }
  return clamp(1.0 - 1.6 * occ, 0.0, 1.0);
}

float rand(vec2 x, float t) {
  x = fract(x * vec2(5.3987, 5.4421));
  x += dot(x.yx, x.xy + vec2(21.5351, 14.3137));
  float xy = x.x * x.y;
  return fract(xy * 95.4307) + fract(xy * 75.04961 + t) - 1.0;
}

float gaussian(float z, float u, float o) {
  return (
    (1.0 / (o * sqrt(TAU))) *
    (exp(-(((z - u) * (z - u)) / (2.0 * (o * o)))))
  );
}

vec3 gaussgrain(float t) {
  vec2 ps = vec2(1.0) / resolution.xy;
  vec2 uv = gl_FragCoord.xy * ps;
  float noise = rand(uv, t);
  noise = gaussian(noise, 0.0, 0.5);
  return vec3(noise);
}

vec4 render(vec2 uv) {
  float camRadius = 1.3;
  vec3 ro = vec3(
    sin(mouselerp.x * PI) * camRadius,
    2.1,
    cos(mouselerp.x * PI) * camRadius);
  vec3 camTarget = vec3(0.0, 0.5, 0.0);
  float camRoll = 0.0;
  mat3 camMatrix = calcLookAtMatrix(ro, camTarget, camRoll);
  vec3 rd = normalize(camMatrix * vec3(uv.x, uv.y, 1.0));
  float d = rayMarch(ro, rd);
  vec3 p = ro + rd * d;
  vec3 sn = normal(p);
  float fd = mapDistance(p);
  float la = light(p, sn, vec3(-3.0, 3.5, 0.0));
  float lb = light(p, sn, vec3(+3.0, 3.5, 0.0));
  float fog = 1.0 / (1.0 + d * d * 0.01 + fd * 5.0);
  vec3 diff = vec3(1.0, 1.0, 1.0) * 0.1;
  float shadowA = getSoftShadow(p, normalize(vec3(-1.0, 2.2, 0.5)), 0.7, 10.0, 1.0);
  float shadowB = getSoftShadow(p, normalize(vec3(+1.0, 2.2, 0.5)), 0.7, 10.0, 1.0);
  float ao = getAmbientOcclusion(p, sn);
  float dp = max(dot(rd, -sn),0.0);
  vec3 col = diff * dp;
  col += la * vec3(1.0, 0.3, 0.5);
  col += lb * vec3(0.2, 0.5, 1.0);
  col = mix(col, col * fog * shadowA * shadowB - gaussgrain(time) * 0.03, 0.75);
  vec3 color = mix(col * ao, vec3(ao), 0.05);
  return vec4(color * 1.2 + vec3(ao) * 0.05, 1.0);
}

vec4 superSample(vec2 uv, float res, int steps, float offset) {
  vec4 rv = vec4(0.0);
  vec2 stepSize = vec2(offset) / res / float(steps);
  uv -= vec2(0.5) / res;
  for (int x = 0; x < steps; ++x) {
    for (int y = 0; y < steps; ++y) {
      vec2 off = vec2(float(x), float(y)) + vec2(0.5);
      rv += render(uv + off * stepSize);
    }
  }
  return rv / float(steps * steps);
}

void main(void) {
  vec2 uv = (gl_FragCoord.xy / resolution.xy * 2.0 - 1.0) * vec2(resolution.x / resolution.y, 1.0);
  float minRes = min(resolution.x, resolution.y);
  #if supersample == 1
    vec4 color = superSample(uv, minRes, 2, 1.25);
  #else
    vec4 color = render(uv);
  #endif
  fragColor = color;
}