Hilbert CA [shadertoy] (basis)

hilbertCA.shadertoy.glsl

// http://bit.ly/supersinfulsilicon

// hilbert CA by Jason Tully is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
// Based on a work at https://www.shadertoy.com/view/ssVXWd
// Permissions beyond the scope of this license may be available at http://bit.ly/supersinfulsilicon

// http://bit.ly/supersinfulsilicon

// hilbert CA by Jason Tully is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
// Based on a work at https://www.shadertoy.com/view/ssVXWd
// Permissions beyond the scope of this license may be available at http://bit.ly/supersinfulsilicon

#define PI (3.141592654f)

float bellcurve(float x) // 0..1 input to 0..1 output
{
    // mu is 0.0 (centered), sigma is 0.5
    const float c = float(PI / -1.25331414f); // optimized magic value - bellcurve perfect match (to six digits of precision) - removes sqrt
                                              // https://www.desmos.com/calculator/xxwdiqa4sk
    x = 2.0f * (2.0f * x - 1.0f);  // converts input range

    return(exp(x * x * c));
}

uint width()
{
    uint mini = uint(min(iResolution.x, iResolution.y));
    
    --mini;
    mini |= (mini >> 1u);
    mini |= (mini >> 2u);
    mini |= (mini >> 4u);
    mini |= (mini >> 8u);
    mini |= (mini >> 16u);
    ++mini;
    
    return(max(1u, mini));
}

float hilbert1D(in uvec2 Position, in uint w){   
    
    uint Index = 0u;
    
    for(uint CurLevel = (w>>1u); CurLevel > 0u; CurLevel >>= 1u){
        uvec2 Region = uvec2(greaterThan((Position & uvec2(CurLevel)), uvec2(0U)));
        Index += CurLevel*CurLevel*((3u*Region.x)^Region.y);
        if(Region.y == 0u){
            if(Region.x == 1u){
                Position = uvec2(w-1u) - Position;
            }
            Position.xy=Position.yx;
        }
    }
    return (float(Index)/float(w*w));
}

vec2 hilbert2D(in uint Index, in uint w){   
    
    uvec2 Position = uvec2(0u);
    
    for(uint CurLevel = 1u; CurLevel < w; CurLevel <<= 1){
        uvec2 Region;
        Region.x = 1u & (Index >> 1u);
        Region.y = 1u & (Index ^ Region.x);
        if(Region.y == 0u){
            if(Region.x == 1u){
                Position = uvec2(CurLevel-1u) - Position;
            }
            Position.xy = Position.yx;
        }
        Position += CurLevel * Region;
        Index >>= 2u;
    }
    return vec2(Position);
}


void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from 0 to 1)
    vec2 uv = fragCoord/iResolution.xy;

    vec4 hn;
    hn.x = textureLod(iChannel0, uv, 0.0f).r;
    hn.y = textureLod(iChannel0, uv, 1.0f).r;
    hn.z = textureLod(iChannel0, uv, 2.0f).r;
    hn.w = textureLod(iChannel0, uv, 4.0f).r;
    
    vec4 nn;
    nn.x = textureLodOffset(iChannel0, uv, 0.0f, ivec2( 0, 1)).r;
    nn.y = textureLodOffset(iChannel0, uv, 0.0f, ivec2( 0,-1)).r;
    nn.z = textureLodOffset(iChannel0, uv, 0.0f, ivec2( 1, 0)).r;
    nn.w = textureLodOffset(iChannel0, uv, 0.0f, ivec2(-1, 0)).r;
    
    float s = 0.0f;
    s = nn.x + nn.y + nn.z + nn.w;
    s = bellcurve(s * 0.25f);
    
    uint w = width();
    
    vec2 hx = hilbert2D(uint(nn.x * float(iFrame)), w);
    vec2 hy = hilbert2D(uint(nn.y * float(iFrame)), w);
    vec2 hz = hilbert2D(uint(nn.z * float(iFrame)), w);
    vec2 hw = hilbert2D(uint(nn.w * float(iFrame)), w);
    
    float h = hilbert1D(uvec2(fragCoord.xy + hx / s + hy / s + hz / s + hw / s), w);

    h = bellcurve(h);
    h = smoothstep(0.9999f, 1.0f, h);
  
    h = h + hn.x * 0.24999f + hn.y * 0.2499f + hn.z * 0.249f + hn.w * 0.248f;
    
    fragColor = vec4(h, hn.xyz);
}