precision highp float;
uniform vec2 resolution;
uniform vec2 mouse;
uniform float time;
uniform sampler2D backbuffer;
out vec4 outColor;

const int maxDepth = 3; // レイの経路の深さの最大値（反射回数+1の値）
const int numSamples = 20; // パストレーシングのサンプル数

const float PI = acos(-1.); // 円周率
const float PI2 = acos(-1.) * 2.;
const float EPS = 0.0001; // レイトレースなどに使う微小量
const float FOV = 60.; // 視野角（degree） 範囲：(0., 180.)
const float fallInterval = 100.; // Voxelが落ちてくる間隔
const float fallSpeed = 4.; // Voxelが落ちる速さ
const float heightRange = 5.; // Voxelの高さの範囲
//const float bumpFactor = 0.7;
const float bumpFactor = 2.;

vec3 camPos; // カメラの座標
vec3 voxelID; // VoxelのID
vec3 voxelNormal; // Voxelの法線ベクトル ((±1, 0, 0) or (0, ±1, 0) or (0, 0, ±1))
vec3 voxelPos; // Voxelの整数座標
float pathSeed = 0.; // パストレーシングで使う乱数のシード

// 1Dの乱数
float hash11(float x) {
  return fract(sin(x) * 43758.5453);
}

// 2Dの乱数
float hash12(vec2 p) {
  return fract(sin(dot(p, vec2(12.9898, 78.233))) * 43758.5453);
}

// 3Dの乱数
float hash13(vec3 p) {
  return fract(sin(dot(p, vec3(12.9898, 78.233, 52.178))) * 43758.5453);
}

// 1Dの乱数(シードを更新)
float random() {
  return hash11(pathSeed++);
}

// 2Dの回転行列
mat2 rotate2D(float a) {
  float s = sin(a), c = cos(a);
  return mat2(c, s, -s, c);
}

// HSVをRGBに変換
vec3 hsv(float h, float s, float v) {
  vec3 res = fract(h + vec3(0, 2, 1) / 3.);
  res = clamp(abs(res * 6. - 3.) - 1., 0., 1.);
  res = (res - 1.) * s + 1.;
  res *= v;
  return res;
}

// 座標pにVoxelが存在するか否か
float map(vec3 p) {
  float res = 1.;
  vec3 ID = floor(p);
  float a = abs(p.x);
  ID.y -= floor(max(a * a * 0.1, 5.) - hash12(ID.xz) * heightRange);
  float T = time * fallSpeed / fallInterval + hash12(ID.xz * 1.852813);
  if(ID.y > 0.){
    float tmp = ID.y / fallInterval + T;
    if(fract(tmp) * fallInterval > 1.) {
      res = 0.;
    }
    ID.y = floor(tmp);
  } else {
    ID.y += floor(T);
  }
  
  p += 0.5 - camPos;
  if(dot(p, p) < 1.) {
    res = 0.;
  }
  voxelID = ID;
  return res;
}

// ref: "Voxel Edges" by iq
// https://www.shadertoy.com/view/4dfGzs
// roから物体表面まで飛ばしたレイの長さ
float castRay(vec3 ro, vec3 rd, const int itr) {
  //vec3 pos = floor(ro);
  voxelPos = floor(ro);
  vec3 ri = 1.0 / rd;
  vec3 rs = sign(rd);
  vec3 dis = (voxelPos - ro + 0.5 + rs * 0.5) * ri; // roからセルの境界までのレイの長さ（X,Y,Z方向）
  
  float res = -1.0;
  vec3 mm = vec3(0);
  vec3 tmp = vec3(0);
  for(int i = 0; i < itr; i++) {
    if(map(voxelPos) > 0.5) {
      res = 1.0;
      break;
    }
    mm = step(dis.xyz, dis.yzx) * step(dis.xyz, dis.zxy); // X,Y,Zのうち、どの方向の隣接セルに進むか
    tmp = mm * rs;
    dis += tmp * ri;
    voxelPos += tmp;
  }
  
  voxelNormal = -tmp;
  vec3 mini = (voxelPos - ro + 0.5 - 0.5 * rs) * ri;
  float t = max(mini.x, max(mini.y, mini.z)); // roからVoxel表面まで伸ばしたレイの長さ
  
  return t * res;
}

// エミッションの色
float emission() {
  return step(hash13(voxelID * 1.552897), 0.08) * 10.;
}

// 物体の色
vec3 objColor(){
  return hsv(hash13(voxelID * 1.347385), 0.6, 1.);
}

// Voxelの凹凸
float bumpFunc(vec3 p, vec3 n) {
  n = abs(n);
  vec2 uv = n.x > 0.5 ? p.yz : n.y > 0.5 ? p.xz : p.xy;
  vec2 ID = floor(uv);
  uv = fract(uv) - 0.5;
  
  if(hash12(ID) < 0.5) {
    uv.y = -uv.y;
  }
  if(uv.y < -uv.x) {
    uv = -uv.yx;
  }
  float d = abs(length(uv - 0.5) - 0.5);
  const float w = 0.15;
  float tmp = w * w - d * d;
  float res = tmp > 0. ? sqrt(tmp) : 0.;
  
  return -res;
}

// ref: "Maze Lattice" by Shane
// https://www.shadertoy.com/view/llGGzh
// Voxelの凹凸ベクトル
vec3 bumpMap(vec3 p, vec3 n) {
  const vec2 e = vec2(.002, 0);
  float ref = bumpFunc(p, n);
  vec3 grad = (vec3(bumpFunc(p - e.xyy, n),
                    bumpFunc(p - e.yxy, n),
                    bumpFunc(p - e.yyx, n)) - ref) / e.x;
  grad -= n * dot(n, grad);
  return normalize(n + grad * bumpFactor);
}

// vを天頂として極座標phi, thetaだけ回転させたベクトル
// ※引数vは長さ1のベクトルである必要あり
vec3 jitter(vec3 v, float phi, float sinTheta, float cosTheta) {
  vec3 xAxis = normalize(cross(v.yzx, v));
  vec3 yAxis = cross(v, xAxis);
  vec3 zAxis = v;
  return (xAxis * cos(phi) + yAxis * sin(phi)) * sinTheta + zAxis * cosTheta;
}

// ref: "パストレーシング - Computer Graphics - memoRANDOM" by Shocker_0x15
// (Japanese article)
// https://rayspace.xyz/CG/contents/path_tracing/

// ref: "GLSL smallpt" by Zavie
// https://www.shadertoy.com/view/4sfGDB

// ref: "[SESSIONS] Syobon's Lobby" by Kamoshika (myself)
// https://www.shadertoy.com/view/ctt3zX

// パストレーシングで得られる色
vec3 pathTrace(vec3 ro, vec3 rd) {
  vec3 acc = vec3(0);
  vec3 mask = vec3(1);
  
  // 最初にカメラからレイを飛ばす
  float t = castRay(ro, rd, 100);
  if(t < 0.) { // レイが物体に衝突しなかった
    return vec3(0.);
  }
  ro += t * rd; // レイの原点を物体表面まで進める
  
  vec3 f = ro - voxelPos - 0.5;
  vec3 n = normalize(f * pow(f, vec3(8.)));
  vec3 bump = bumpMap(ro, voxelNormal);
  n = normalize(n + bump);
  
  // ※本来はサンプル数の分だけカメラからレイを飛ばすため、numSamplesをかける
  acc += mask * emission() * float(numSamples);
  mask *= objColor();
  
  // 次に、物体表面からランダムにレイを飛ばす
  vec3 ro0 = ro + n * EPS;
  vec3 n0 = n;
  vec3 mask0 = mask;
  for(int i = 0; i < numSamples; i++) {
    ro = ro0;
    n = n0;
    mask = mask0;
    for(int depth = 1; depth < maxDepth; depth++) {
      float ur = random(); // 一様乱数
      // 重点的サンプリングを行うため、半球面内cos分布を使用する
      rd = jitter(n, random() * PI2, sqrt(1. - ur), sqrt(ur)); // 次のレイの方向
            
      t = castRay(ro, rd, 15);
      if(t < 0.) { // レイが物体に衝突しなかった
        break;
      }
      ro += t * rd;
      
      f = ro - voxelPos - 0.5;
      n = normalize(f * pow(f, vec3(8.)));
      bump = bumpMap(ro, voxelNormal);
      n = normalize(n + bump);
      
      acc += mask * emission();
      mask *= objColor();
      ro += n * EPS; // 現在の物体表面を避けるために、少し浮かせる
    }
  }
  
  acc /= float(numSamples);
  acc = clamp(acc, 0., 1.);
  
  return acc;
}

void main() {
  vec2 uv = (gl_FragCoord.xy * 2. - resolution) / min(resolution.x, resolution.y); // 座標の正規化
  
  camPos = vec3(0.5, 7.5, -time * 2.); // カメラの位置（レイの原点）
  vec3 dir = normalize(vec3(0.5, -0.3, -1)); // カメラの向き
  dir.xy *= rotate2D(time * 0.2);
  vec3 side = normalize(cross(dir, vec3(0, 1, 0)));
  vec3 up = cross(side, dir);
  vec3 rd = normalize(uv.x * side + uv.y * up + dir / tan(FOV / 360. * PI)); // レイの向き
  
  // パストレーシングで使う乱数のシードを初期化する
  float T = fract(time / 10.) * 500.;
  pathSeed = hash12(gl_FragCoord.xy * PI) * 500.;
  pathSeed += hash12(gl_FragCoord.xy + pathSeed + T) * 500.;
  
  // パストレーシングをする
  vec3 col = vec3(0);
  col += pathTrace(camPos, rd);
  col = pow(col, vec3(1. / 2.2)); // ガンマ補正
  
  // パストレーシングで得られた色の分散を低減するために前のフレームの色を合成する
  float tex = 0.6;
  col = mix(col, texture(backbuffer, gl_FragCoord.xy / resolution).rgb, tex);
  
  outColor = vec4(col, 1.);
}