#define PI 3.141593
#define MAX_VOXELS 100
#define SECONDARY_STEPS 10
#define SAMPLES_PER_FRAME 30
#define BOUNCES 3
#define FDIST 0.5
#define ATMOSPHERE_TRANSMISSION 0.99
#define SCATTER_FACTOR 0.5
#define EPS 0.01

#define PERIOD 11.
#define BUILDINGRAD 4.

struct Hit {
    float t;
    int mat;
    vec3 n;
    vec3 id;
};

float noise2d(in vec2 uv) {
    return fract(814.*sin(uv.x*15829.+uv.y*874.));
}

uint seed;
uint wang_hash()
{
    seed = uint(seed ^ uint(61)) ^ uint(seed >> uint(16));
    seed *= uint(9);
    seed = seed ^ (seed >> 4);
    seed *= uint(0x27d4eb2d);
    seed = seed ^ (seed >> 15);
    return seed;
}

float GetRandom()
{
    return float(wang_hash()) / 4294967296.0;
}

vec2 GetRandom2(){return vec2(GetRandom(), GetRandom());}

vec3 randnorm() {
    vec2 utheta = GetRandom2() * 2. - 1.;
    utheta.y *= PI;
    float rho = sqrt(1.-utheta.x * utheta.x);
    return vec3(cos(utheta.y)*rho, sin(utheta.y)*rho, utheta.x);
}


vec2 condmin(in vec2 d1, in vec2 d2) {
    return vec2(min(d1.x, d2.x), mix(d1.y, d2.y, step(d2.x, d1.x)));
}

//occupancy function + material
float occupancy(in vec3 id, out int mat) {
    vec2 block = floor(id.xy/PERIOD+0.5);
    vec3 modid = vec3(mod(id.xy+0.5*PERIOD, PERIOD)-0.5*PERIOD, id.z);
    float randval1 = noise2d(block);
    float randval2 = noise2d(vec2(block.y, randval1));

    vec2 mindist = vec2(max(modid.z-1., PERIOD*0.5 - length(modid.xy)), 3.); //base
    vec2 buildings = vec2(max(length(modid.xy)-BUILDINGRAD, -30.+modid.z-100.*randval2), 2.);
    buildings.x = max(buildings.x, -length(vec2(mix(modid.x, modid.y, mod(floor(modid.z/4.), 2.0)), mod(modid.z, 4.)-2.)) + 2.);
    mindist = condmin(mindist, buildings);
    mindist = mix(mindist, vec2(modid.z, 1.), step(0.3, randval1));
    float occ = step(mindist.x, 0.5);
    mat = int(occ * round(mindist.y));
    return occ;
}

Hit voxtrace(in vec3 ro, in vec3 rd, int iters) {
    Hit h;
    h.mat = 0;
    h.t = 0.;
    // box marching with fb39ca4's DDA
    h.id = floor(ro);
    vec3 ri = 1.0/rd;
    vec3 rs = sign(rd);
    vec3 dis = (h.id-ro + 0.5 + rs*0.5) * ri;
    vec3 mm = vec3(0.0);
    for (int i=0; i<iters; i++) {
        mm = step(dis.xyz, dis.yzx) * step(dis.xyz, dis.zxy);
		dis += mm * rs * ri;
        h.id += mm * rs;
        if (occupancy(h.id, h.mat) > 0.5) {
            break;
    	}
    }

    h.n = -mm*rs;

	vec3 mini = (h.id-ro + 0.5 - 0.5*vec3(rs))*ri;
	h.t = max ( mini.x, max ( mini.y, mini.z ) );

    return h;
}

vec3 bgcol(in vec3 rd) {
    return mix(vec3(0., 0., 1.), vec3(0.6, 0.8, 1.), 1.-pow(abs(rd.z), 2.));
}

vec3 pathtrace(in vec3 eye, in vec3 rd) {
    vec3 col = vec3(0.);
    vec3 thru = vec3(1.);
    for (int i=0; i<BOUNCES; i++) {
        Hit h = voxtrace(eye, rd, i == 0 ? MAX_VOXELS : SECONDARY_STEPS);
        //thru *= pow(ATMOSPHERE_TRANSMISSION, h.t);
        if (h.mat == 0) {
            col += thru * bgcol(rd);
            break;
        }
        eye = eye + h.t * rd;
        vec3 nr = randnorm();
        vec3 nref = reflect(rd, h.n);
		nr = mix(nref, normalize(h.n + nr), SCATTER_FACTOR);
        vec3 albedo = vec3(1.);
        vec3 emissive = vec3(0.);
        if (h.mat == 1) {
            albedo = vec3(0.3, 0.7, 0.9);
        } else if (h.mat == 2) {
            albedo = vec3(0.6, 0.9, 0.5);
        } else if (h.mat == 3) {
            albedo = vec3(1.);
            //emissive = vec3(0.5, 0.4, 0.2);
        }

        col += thru * emissive;
        thru *= albedo;

        float maxthru = max(thru.x, max(thru.y, thru.z));
        if (GetRandom() > maxthru) {
            break;
        }
        thru *= 1.0/maxthru;

        rd = nr;
        eye += EPS * rd;
    }
    return col;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    seed = uint(uint(fragCoord.x) * uint(1973) + uint(fragCoord.y) * uint(9277) + uint(iFrame) * uint(26699)) | uint(1);
    float mouseY = iMouse.y < 1. ? 0. : (0.5-iMouse.y/iResolution.y) * PI;
    float mouseX = iMouse.x < 1. ? iTime*0.25 : -(iMouse.x/iResolution.x) * 2. * PI;
    vec2 uv = (fragCoord - 0.5*iResolution.xy)/iResolution.x;
	vec3 eye = vec3(PERIOD*0.5, iTime * 9., 6.1);
    vec3 w = vec3(cos(mouseX) * cos(mouseY), sin(mouseX) * cos(mouseY), -sin(mouseY));
    vec3 u = normalize(cross(w, vec3(0., 0., 1.)));
    vec3 v = cross(u, w);
    vec3 rd = normalize(FDIST*w + uv.x*u + uv.y*v);

    vec3 col = vec3(0.);
    float weight = 0.0;
    for (int i=0; i<SAMPLES_PER_FRAME; i++) {
        float newweight = weight + 1.0;
        vec3 newcol = pathtrace(eye, rd);
        col = (col * weight + newcol)/newweight;
        weight = newweight;
    }

    //dummy visualization
    fragColor = vec4(pow(col, vec3(0.45)), 1.0);
    //fragColor = vec4((sin(h.t*10.)*.1+.5) * float(h.mat) * (h.n*.5 +.5), 1.);
    //fragColor = vec4(fract(h.id/10.), 1.);
    //fragColor = vec4(vec3(h.mat)/3., 1.);
}