// with small optimal tweakings from @nothings

// WIP
//#define CLOSED_FORM

const float tau = 6.283185;

bool cube(vec3 org, vec3 dir, out float near, out float far) {
    vec3 p = -org/dir;
    vec3 q = 0.5/abs(dir);
	vec3 tmin = p - q;
	vec3 tmax = p + q;
    near = max(tmin.x,max(tmin.y,tmin.z));
	far = min(tmax.x,min(tmax.y,tmax.z));
	return near < far && far > 0.0;
}


void transformray (vec3 ro, vec3 rd, mat2 rotationY, vec3 offset, out vec3 outro, out vec3 outrd)
{
	outro = ro + offset;
	outro = vec3(rotationY * outro.xz, outro.y).xzy;
	outrd = vec3(rotationY * rd.xz, rd.y).xzy;
}

void trilinear_integral_factors(out float f[8], vec3 ro, vec3 rd, float t) {
    float t2 = 0.5*t*t;
    float t3 = (2.0/3.0)*t2*t;
    float t4 = t2*t2;

    float dxdydz = rd.x*rd.y*rd.z*t4;
    float ddo = (rd.y*rd.z*ro.x + rd.x*rd.z*ro.y + rd.x*rd.y*ro.z)*t3;
    float doo = (rd.z*ro.x*ro.y + rd.y*ro.x*ro.z + rd.x*ro.y*ro.z)*t2;
    float oxoyoz = ro.x*ro.y*ro.z*t;
    float dxdydz_ddo_doo_oxoyoz = dxdydz + ddo + doo + oxoyoz;

    float dxdy = rd.x*rd.y*t3;
    float dxdz = rd.x*rd.z*t3;
    float dydz = rd.y*rd.z*t3;

    float oxdy = rd.y*ro.x*t2;
    float oxdz = rd.z*ro.x*t2;

    float dxoy = rd.x*ro.y*t2;
    float oydz = rd.z*ro.y*t2;

    float dxoz = rd.x*ro.z*t2;
    float dyoz = rd.y*ro.z*t2;

    float dx = rd.x*t2;
    float dy = rd.y*t2;
    float dz = rd.z*t2;

    float oxoy = ro.x*ro.y*t;
    float oxoz = ro.x*ro.z*t;
    float oyoz = ro.y*ro.z*t;

    float ox = ro.x*t;
    float oy = ro.y*t;
    float oz = ro.z*t;

    f[0] = ( // 000
    - dxdydz_ddo_doo_oxoyoz
    + dxdy + dxdz + dydz
    + oxdy + oxdz + dxoy + oydz + dxoz + dyoz - dx - dy - dz
    + oxoy + oxoz + oyoz - ox - oy - oz + t
    );
    f[1] = ( // 001
    + dxdydz_ddo_doo_oxoyoz
    - dxdy - dxdz
    - oxdy - oxdz - dxoy - dxoz + dx
    - oxoy - oxoz + ox
    );
    f[2] = ( // 010
    + dxdydz_ddo_doo_oxoyoz
    - dxdy - dydz
    - oxdy - dxoy - oydz - dyoz + dy
    - oxoy - oyoz + oy
    );
    f[4] = ( // 100
    + dxdydz_ddo_doo_oxoyoz
    - dxdz - dydz
    - oxdz - oydz - dxoz - dyoz + dz
    - oxoz - oyoz + oz
    );
    f[6] = ( // 110
    - dxdydz_ddo_doo_oxoyoz
    + dydz
    + oydz + dyoz
    + oyoz
    );
    f[5] = ( // 101
    - dxdydz_ddo_doo_oxoyoz
    + dxdz
    + oxdz + dxoz
    + oxoz
    );
    f[3] = ( // 011
    - dxdydz_ddo_doo_oxoyoz
    + dxdy
    + oxdy + dxoy
    + oxoy
    );
    f[7] = ( // 111
    + dxdydz_ddo_doo_oxoyoz
    );

}

float dot8(float c[8], float f[8]) {
    return
          c[0] * f[0]
        + c[1] * f[1]
        + c[2] * f[2]
        + c[3] * f[3]
        + c[4] * f[4]
        + c[5] * f[5]
        + c[6] * f[6]
        + c[7] * f[7];
}

const float A = 1.0;
const float B = 0.0;
vec4 C[8] = vec4[](
#if 1
    vec4(0,0,0,B),
    vec4(1,0,0,A),
    vec4(0,1,0,A),
    vec4(1,1,0,B),
    vec4(0,0,1,A),
    vec4(1,0,1,B),
    vec4(0,1,1,B),
    vec4(1,1,1,A)
#elif 0
    vec4(1,1,1,A),
    vec4(1,1,1,A),
    vec4(1,1,1,A),
    vec4(1,1,1,A),
    vec4(1,1,1,B),
    vec4(1,1,1,B),
    vec4(1,1,1,B),
    vec4(1,1,1,B)
#elif 1
    vec4(0,0,0,B),
    vec4(1,0,0,B),
    vec4(0,1,0,B),
    vec4(1,1,0,B),
    vec4(0,0,1,A),
    vec4(1,0,1,A),
    vec4(0,1,1,A),
    vec4(1,1,1,A)
#else
    vec4(0,0,0,B),
    vec4(1,0,0,A),
    vec4(0,1,0,A),
    vec4(1,1,0,A),
    vec4(0,0,1,A),
    vec4(1,0,1,A),
    vec4(0,1,1,A),
    vec4(1,1,1,B)
#endif
);

// smallest non-inf number: exp(-103.9) or exp2(-149.99999)

vec4 tf(vec4 s) {
    return s;
}

vec4 tfinv(vec4 s) {
    return s;
}

vec4 map(vec3 p) {
    vec4 c00 = mix(tf(C[0]), tf(C[1]), p.x);
    vec4 c01 = mix(tf(C[2]), tf(C[3]), p.x);
    vec4 c10 = mix(tf(C[4]), tf(C[5]), p.x);
    vec4 c11 = mix(tf(C[6]), tf(C[7]), p.x);
    vec4 c0 = mix(c00, c01, p.y);
    vec4 c1 = mix(c10, c11, p.y);
    vec4 c = mix(c0, c1, p.z);
    return c;
}

// 8 is the lowest that still looks okay
// there's only a faint visible difference between 32 and 256
// particularly the peaks become more intense at low iterations
const int N = 256;
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
	// camera
	vec2 q = fragCoord.xy/iResolution.xy;
	vec2 p = -1.0 + 2.0 * q;
	p.x *= iResolution.x/iResolution.y;
    p *= 0.55; // move camera closer
	vec3 camro = normalize(vec3(1.0, -0.1 + sin(iTime)*0.5+0.5, 0.0));
	vec3 w = -camro;
	camro *= 2.5;
	vec3 u = normalize(cross( vec3(0.0, 1.0, 0.0), w ));
	vec3 v = normalize(cross(w,u));
	vec3 camrd = normalize(p.x * u + p.y * v + 1.5 * w);
	fragColor = vec4(0.0);

	// rotation
	float angle = 5.0 * iMouse.x / iResolution.x;
	if( iMouse.z < 0.5 )
		angle = iTime + 4.7;
	float ca = cos(angle);
	float sa = sin(angle);
	mat2  m = mat2(ca, -sa, sa, ca);

	float far, near;
	vec3 ro, rd;

    vec3 c = vec3(0.0);
#ifdef CLOSED_FORM
    float k = 1.0;
#else
    float k = 1.0;
#endif
	// cube
	transformray(camro, camrd, m, vec3(0.0), ro, rd );
    if (cube (ro, rd, near, far)) {
        float d = (far - near);
        vec3 p0 = (ro + rd * near)+0.5;
        vec3 p1 = (ro + rd * far)+0.5;
#ifdef CLOSED_FORM
        float f[8];
        vec3 rd = p1 - p0;
        float l = length(rd);
        trilinear_integral_factors(f, p0, normalize(rd), l);
        //trilinear_integral_factors(f, p0, rd, 1.0);
        k = 0.0;
        float CF = 0.99;
        for (int i = 0; i < 8; ++i) {
            k += log(1.0 - C[i].w*CF)*f[i];
            c += log(1.0 - C[i].rgb*C[i].w*CF)*f[i];
        }
        k = exp(k);
        c = 1.0 - exp(c);
#else
        float s = 1.0/float(N);
        vec4 c0 = map(p0);
        vec4 c1 = map(p1);
        // solve front to back
        for (int i = 0; i <= N; ++i) {
            float x = float(i)*s;
//#define CLOSED
#ifdef CLOSED
        float a = c0.w;
        float b = c1.w;
        //float k0 = exp(-D*log(a)/(a - b));
        float k0 = d*x;
        float k1 = log(a + (b-a)*x)*d/(a - b);
        float k1b = log(a + (b-a)*x)*(d/(a - b) - 1.0);
        float k2 = -log(a)*d/(a-b);
        vec3 u = c0.rgb;
        vec3 v = c1.rgb;
#endif
            vec3 p = clamp(p0 + (p1 - p0)*x, vec3(0.0), vec3(1.0));
     		vec4 src = tfinv(map(p));
            vec3 col = src.rgb; // color
            // two possible volume curves, experimentally found
            //float tl = 1.0/(1.0 + exp(1.0/(1.0 - a) - 1.0/a)); // planck-taper function

            // front to back alpha blending iteration
#ifdef CLOSED
            float C = b-a;
            vec3 D = v-u;
            float Tl = a + C*x;
            k = clamp(exp(d*x)*pow(1.0-x*C/(1.0-a), d/C),0.0,1.0);
            //c = c + s*d*(u + x*D)*(a + C*x)*k/(1.0 - (a + C*x));

            vec3 T0 = u*a;
            vec3 T1 = (D*a + u*C)*x;
            vec3 T2 = D*C*x*x;

            c = c + s*d*(T0 + T1 + T2)*k/(1.0 - (a + C*x));

#else
            float a = src.w; // normalized loglevel
            float tl = 1.0 - 1.0/(1.0 - a);
            float ttl = pow(exp(tl), d*s); // distribute over s=1/N iterations
            c = c + col*(1.0 - ttl)*k;
			k *= ttl; // compute remaining visibility
#endif
    	}
#endif
    }
#ifdef CLOSED_FORM
    c = c + (texture(iChannel0, p*2.0).rrr * 0.2 + 0.2)*k;
#else
    c = c + (texture(iChannel0, p*2.0).rrr * 0.2 + 0.2)*k;
    //c = vec3(k);
#endif
#if 0
    if (max(c.r,max(c.g,c.b)) > 1.0) {
        c = vec3(1.0,0.0,0.0);
    }
#endif
    fragColor = vec4(pow(c, vec3(1.0/2.2)), 1.0);
}