#define MDIST 150.0
#define STEPS 300.0
#define pi 3.1415926535
#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))
#define fs(i) (fract(sin((i)*114.514)*1919.810))
//iq palette
vec3 pal( in float t, in vec3 a, in vec3 b, in vec3 c, in vec3 d ){
    return a + b*cos(2.*pi*(c*t+d));
}
float h21 (vec2 a) {
    return fract(sin(dot(a.xy,vec2(12.9898,78.233)))*43758.5453123);
}
float h11 (float a) {
    return fract(sin((a)*12.9898)*43758.5453123);
}
float box(vec3 p, vec3 b){
    vec3 d = abs(p)-b;
    return max(d.x,max(d.y,d.z));
}
float ebox(vec3 p, vec3 b){
  vec3 q = abs(p) - b;
  return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}

//Based on code from bigwings comment here
//https://www.shadertoy.com/view/Wl3fD2
float dibox(vec3 p,vec3 b,vec3 rd){
    p/=b;
    vec3 dir = sign(rd)*.5;
    vec3 rc = (dir-p)/rd;
    rc*=b;
    float dc = min(min(rc.x, rc.y), rc.z)+0.001;
    return dc;
}

vec3 rdg = vec3(0);

struct RectSubdivResult {
    float volume;
    vec3 center;
    vec3 dimension;
    float id;
};

/**
 * Ref: https://www.shadertoy.com/view/7sKGRy
 * @param p The input position
 * @param scale scale of the domain of the fractal
 */
RectSubdivResult rectSubdiv( vec3 p, vec3 scale ) {
    float t = iTime;

    // several constants
    const int ITERS = 12;
    const int MIN_ITERS = 1;
    const float MIN_SIZE = 0.15;
    const float BREAK_CHANCE = 0.0;
    const float PAD_FACTOR = 1.01;

    // the domain of the fractal being generated
    // will be modified in the iteration part
    vec3 domainMin = vec3( -0.5 ) * scale;
    vec3 domainMax = vec3( 0.5 ) * scale;

    // id of the individual cube in the fractal
    float id = 0.0;

    // random seed of cut positions
    float seed = floor( t / 6.0 ) + 0.1;

    // size of the current box determined by domainMin / domainMax
    vec3 dimension = domainMax - domainMin;

    for ( int i = 0; i < ITERS; i ++ ) {
        float fi = float( i );

        // divide the box into eight
        vec3 divideHash = vec3(
            h21( vec2( fi + id, seed ) ),
            h21( vec2( fi + id + 2.44, seed ) ),
            h21( vec2( fi + id + 7.83, seed ) )
        );
        vec3 divide = divideHash * dimension + domainMin;

        // let the division line cut the box not too thin
        divide = clamp( divide, domainMin + MIN_SIZE * PAD_FACTOR, domainMax - MIN_SIZE * PAD_FACTOR );

        // does this cut the box to the minimum preferrable size?
        vec3 minSizeOfAxis = min( abs( domainMin - divide ), abs( domainMax - divide ) );
        float minSize = min( minSizeOfAxis.x, min( minSizeOfAxis.y, minSizeOfAxis.z ) );
        bool isSmallEnough = minSize < MIN_SIZE;

        bool willBreak = false;
        if ( i - 1 > MIN_ITERS && h11( id ) < BREAK_CHANCE ) { willBreak = true; }
        if ( isSmallEnough && i - 1 > MIN_ITERS || i == ITERS - 1 ) { willBreak = true; }
        if( willBreak ) {
            // id = i * 0.1 * seed;
            break;
        }

        // update the box domain
        domainMax = mix( domainMax, divide, step( p, divide ) );
        domainMin = mix( divide, domainMin, step( p, divide ) );

        // id will be used for coloring and hash seeding
        vec3 diff = mix( -divide, divide, step( p, divide ) );
        id = length( diff + 10.0 );

        // recalculate the dimension
        dimension = domainMax - domainMin;
    }

    // calculate volume and center of the box
    float volume = dimension.x * dimension.y * dimension.z;
    vec3 center = ( domainMin + domainMax ) / 2.0;

    // prepare the result
    RectSubdivResult result;
    result.volume = volume;
    result.center = center;
    result.dimension = dimension;
    result.id = id;

    return result;
}

vec2 blocks( vec3 p, vec3 scale, vec3 rd ) {
    float t = iTime;

    const float MAX_CENTER_DIST = 4.5;
    const float MAX_VOLUME = 5.0;
    const float DESTRUCTION_CHANCE = 0.5;

    RectSubdivResult rectSubdivResult = rectSubdiv( p, scale );

    float volume = rectSubdivResult.volume;
    vec3 center = rectSubdivResult.center;
    vec3 dimension = rectSubdivResult.dimension;
    float id = rectSubdivResult.id;

    //huge improvment in performance by using distance to intersection of empty cell
    //to remove boxes (instead of using a negative box sdf)
    //But it seems to cause artifacts rarely, idk why
    float b = dibox( p - center, dimension, rd );

    float shr = 1.0-abs(pow(abs(cos(mod(t,6.)*pi/6.)),6.0));
    shr = smoothstep(0.,1.,shr);
    vec3 d = abs(center);
    center.y -= dimension.y * ( 1.0 - shr ) * 0.5;
    dimension.y *= shr;
    float a = box( p - center, dimension * 0.5 );

    //I found this helps to remove some of the artifacts from using the empty box intersection
    if( abs( p.x ) > scale.x * 0.5 ) { b = -a; }
    if( abs( p.z ) > scale.z * 0.5 ) { b = -a; }

    a = min( a, b );
    if( max( d.x, max( d.y * 0.5, d.z ) ) > MAX_CENTER_DIST ) { a = b; }
    else if ( volume > MAX_VOLUME ) { a = b; }
    else if ( h11( id * 1.1 ) < DESTRUCTION_CHANCE ) { a = b; }

    id = h11(id)*1000.0;

    return vec2(a,id);
}

vec2 map(vec3 p){
    float t = iTime;
    vec3 po = p;
    vec2 a = vec2(1);
    vec3 scl = vec3(10.0,10.,10);
    vec3 rd2 = rdg;
    a = blocks(p,scl,rdg)+0.02;

    a.x = max(box(p,vec3(scl*0.49)),a.x);

    return a;
}
vec3 norm(vec3 p){
    vec2 e = vec2(0.00005,0.);
    return normalize(map(p).x-vec3(
    map(p-e.xyy).x,
    map(p-e.yxy).x,
    map(p-e.yyx).x));
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord-0.5*iResolution.xy)/iResolution.y;
    float t = iTime;
    vec3 col = vec3(0);
    vec3 ro = vec3(0,5,-20);
    if(iMouse.z>0.){
    ro.zx*=rot(7.0*iMouse.x/iResolution.x);
    }
    else ro.zx*=rot(iTime*0.3);
    vec3 lk = vec3(0,0.,0);
    vec3 f = normalize(lk-ro);
    vec3 r = normalize(cross(vec3(0,1,0),f));
    vec3 rd = normalize(f*0.95+uv.x*r+uv.y*cross(f,r));
    rdg = rd;
    vec3 p = ro;
    float dO = 0.;
    vec2 d = vec2(0);
    bool hit = false;
    for(float i = 0.; i<STEPS; i++){
        p = ro+rd*dO;
        d = map(p);
        dO+=d.x*0.99;
        if(abs(d.x)<0.0001){
            hit = true;
            break;
        }
        if(d.x>MDIST){
            dO=MDIST;
            break;
        }
    }
    if(hit){
        vec3 ld = normalize(vec3(0.5,1,-1));
        vec3 n = norm(p);
        vec3 r = reflect(rd, n);
        vec3 e = vec3(0.5);

        vec3 al = pal(d.y*0.1,e*1.2,e,e*2.0,vec3(0,0.33,0.66));
        if(d.y==2.0) al = vec3(1.);
        col = al;

        float diff = length(sin(n*2.)*.5+.8)/sqrt(3.);
        col = al*diff;

        float shadow = 1.;
        rdg = ld;
        for(float h = 0.05; h<50.;){
            float dd = map(p+ld*h).x;
            if(dd<0.001){shadow = 0.6; break;}
            h+=dd;
        }
        col*=shadow;
    }
    vec3 bg = mix(vec3(0.173,0.231,0.686),vec3(0.361,0.753,1.000),rd.y*0.5+0.5);
    col = mix(col,bg,dO/MDIST);
    fragColor = vec4(col,1.0);
}
/*
#define AA 2.0
#define ZERO min(0.0,iTime)
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    float px = 1.0/AA;
    vec4 col = vec4(0);

    if(AA==1.0) {render(col,fragCoord); fragColor = col; return;}

    for(float i = ZERO; i <AA; i++){
        for(float j = ZERO; j <AA; j++){
            vec4 col2;
            vec2 coord = vec2(fragCoord.x+px*i,fragCoord.y+px*j);
            render(col2,coord);
            col.rgb+=col2.rgb;
            rdg = vec3(0);
        }
    }
    col/=AA*AA;
    fragColor = vec4(col);
}
*/