#define PI 3.14159

#define VOXEL_NONE  0
#define VOXEL_WATER 1
#define VOXEL_SAND  2
#define VOXEL_EARTH 3
#define VOXEL_STONE 4
#define VOXEL_GRASS 5

// from https://www.shadertoy.com/view/4sfGzS
float noise( in vec3 x )
{
    vec3 p = floor(x);
    vec3 f = fract(x);
	f = f*f*(3.0-2.0*f);

    vec2 uv = p.xy + f.xy;
	vec2 rg = vec2(texture( iChannel0, (uv+vec2(37.0,17.0)*p.z+0.5)/256.0, -100.0 ).x,
                   texture( iChannel0, (uv+vec2(37.0,17.0)*(p.z+1.0)+0.5)/256.0, -100.0 ).x );
	return mix( rg.x, rg.y, f.z );
}

float map(vec2 p)
{
    float h = texture(iChannel0, p/800.0).r;
    return h;
}

float mmap(vec2 p)
{
    float f = map(p) * 0.3;
    return texture(iChannel2, vec2(f, 0.1)).r;
}

int getVoxelAt(ivec3 ip)
{
    if (ip.y <= 0) return VOXEL_WATER;

    float h = map(vec2(ip.xz))*20.0-2.0;

    if (float(ip.y) <= h)
        return VOXEL_SAND;

    return VOXEL_NONE;
}

float dfVoxel(vec3 p, int voxelType)
{
    return length(p) - 0.5;
}

vec3 nrmVoxel(vec3 p, int voxelType)
{
    vec2 dd = vec2(0.001,0.0);
    float base = dfVoxel(p, voxelType);
    return normalize(vec3(
        dfVoxel(p+dd.xyy, voxelType) - base,
        dfVoxel(p+dd.yxy, voxelType) - base,
        dfVoxel(p+dd.yyx, voxelType) - base
    ));
}

void voxelMarch(vec3 ro, vec3 rd, out ivec3 hitVoxels[4], out int hitCount, out float mAccums[4])
{
    hitCount = 0;

    ivec3 mapPos = ivec3(floor(ro));
    vec3 deltaDist = abs(vec3(length(rd)) / rd);
    ivec3 rayStep = ivec3(sign(rd));
    vec3 sideDist = (sign(rd) * (vec3(mapPos) - ro) + (sign(rd) * 0.5) + 0.5) * deltaDist;
    bvec3 mask;

    float mAccum = 0.0;

    for (int i = 0; i < 128; i++) {

        // check current position for voxel
        if (getVoxelAt(mapPos) != VOXEL_NONE) {
            // no non-const indexing? :<
            if (hitCount == 0) hitVoxels[0] = mapPos;
            else if (hitCount == 1) hitVoxels[1] = mapPos;
            else if (hitCount == 2) hitVoxels[2] = mapPos;
            else if (hitCount == 3) hitVoxels[3] = mapPos;
            hitCount++;
            if (hitCount == 4) return;
        }

        mAccum += mmap(vec2(mapPos.xz));
        mAccums[hitCount] = mAccum;

        // march forward to next position by discrete digital analyzer
        if (sideDist.x < sideDist.y) {
            if (sideDist.x < sideDist.z) {
                sideDist.x += deltaDist.x;
                mapPos.x += rayStep.x;
                mask = bvec3(true, false, false);
            } else {
                sideDist.z += deltaDist.z;
                mapPos.z += rayStep.z;
                mask = bvec3(false, false, true);
            }
        } else {
            if (sideDist.y < sideDist.z) {
                sideDist.y += deltaDist.y;
                mapPos.y += rayStep.y;
                mask = bvec3(false, true, false);
            } else {
                sideDist.z += deltaDist.z;
                mapPos.z += rayStep.z;
                mask = bvec3(false, false, true);
            }
        }
    }
}

void resolveHitVoxels(
    vec3 ro, vec3 rd, ivec3 hitVoxels[4], float mAccums[4], int hitCount,
    out ivec3 hitVoxel, out vec3 hit, out int terrainType, out float mAccum)
{
  mAccum = mAccums[0];
  for (int i=0; i<4; i++) {
    if (i == hitCount) return;

    hitVoxel = hitVoxels[i];
    mAccum = mAccums[i];
    terrainType = getVoxelAt(hitVoxel);
    vec3 hitVoxelCenter = vec3(hitVoxel) + 0.5;

    // intersect with voxel cube
    vec3 cubeIntersect = (hitVoxelCenter - ro - 0.5*sign(rd))/rd;
    float d = max(cubeIntersect.x, max(cubeIntersect.y, cubeIntersect.z));

    // fallback in case of no distance intersection
    hit = ro + rd * (d - 0.01) - hitVoxelCenter;

    // attempt better intersect with distance marching
    float diff;
    vec3 p = ro + rd * d;
    for (int j=0; j<4; j++) {
      diff = dfVoxel(p - hitVoxelCenter, terrainType);
      d += diff;
      p = ro + rd * d;
    }
    if (diff < 0.05) { // good enough distance marched intersection
      hit = p - hitVoxelCenter;
      return;
    }
  }
}

vec3 doColoring(vec3 hit, int terrainType, vec3 ldir)
{
    vec3 n = nrmVoxel(hit, terrainType);
    float diffuse = max(dot(-ldir, n), 0.1);

    // render
    vec3 color = vec3(0.0);
    if (terrainType == VOXEL_WATER) color = vec3(0.6, 0.6, 1.0);
    if (terrainType == VOXEL_EARTH) color.r = 1.0;
    if (terrainType == VOXEL_SAND) color.rg = vec2(1.0, 0.6);
    if (terrainType == VOXEL_STONE) color.rgb = vec3(0.5);
    if (terrainType == VOXEL_GRASS) color.g = 1.0;

    color *= (0.8 + 0.2*diffuse);

    return color;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // camera stolen from Shane :) https://www.shadertoy.com/view/ll2SRy
	vec2 uv = (fragCoord - iResolution.xy*.5 )/iResolution.y;
    vec3 rd = normalize(vec3(uv, (1.-dot(uv, uv)*.5)*.5));
    vec3 ro = vec3(0., 20., iTime*3.0);
    float t = sin(iTime * 0.2) + noise(ro/32.0);
    ro.y += 4.0*t;
	float cs = cos( t ), si = sin( t );
    rd.yz = mat2(cs, si,-si, cs)*rd.yz;
    rd.xz = mat2(cs, si,-si, cs)*rd.xz;

    // voxel march into the scene storing up to four intersections
    int hitCount;
    ivec3 hitVoxels[4];
    float mAccums[4];
    voxelMarch(ro, rd, hitVoxels, hitCount, mAccums);

    // resolve to one accurate intersection by distance marching
    int terrainType = VOXEL_NONE;
    ivec3 hitVoxel;
    vec3 hit;
    float mAccum;
    resolveHitVoxels(ro, rd, hitVoxels, mAccums, hitCount, hitVoxel, hit, terrainType, mAccum);

    vec3 hitGlobal = vec3(hitVoxel) + hit;

    // color
    vec3 ldir = normalize(hitGlobal - ro);
    vec3 color = doColoring(hit, terrainType, ldir);

    // fog
    float fog = smoothstep(1.0, 0.0, length(hitGlobal - ro)/96.0);
    color *= fog;

    color += mAccum/32.0;

	fragColor = vec4(color,1.0);
}