Walk-on-spheres

infill.glsl

// See https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/

// Random numbers using code at
// https://stackoverflow.com/a/17479300

// A single iteration of Bob Jenkins' One-At-A-Time hashing algorithm.
uint hash( uint x ) {
    x += ( x << 10u );
    x ^= ( x >>  6u );
    x += ( x <<  3u );
    x ^= ( x >> 11u );
    x += ( x << 15u );
    return x;
}


// Compound versions of the hashing algorithm I whipped together.
uint hash( uvec2 v ) { return hash( v.x ^ hash(v.y)                         ); }
uint hash( uvec3 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z)             ); }
uint hash( uvec4 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z) ^ hash(v.w) ); }


// Construct a float with half-open range [0:1] using low 23 bits.
// All zeroes yields 0.0, all ones yields the next smallest representable value below 1.0.
float floatConstruct( uint m ) {
    const uint ieeeMantissa = 0x007FFFFFu; // binary32 mantissa bitmask
    const uint ieeeOne      = 0x3F800000u; // 1.0 in IEEE binary32

    m &= ieeeMantissa;                     // Keep only mantissa bits (fractional part)
    m |= ieeeOne;                          // Add fractional part to 1.0

    float  f = uintBitsToFloat( m );       // Range [1:2]
    return f - 1.0;                        // Range [0:1]
}


// Pseudo-random value in half-open range [0:1].
float random( float x ) { return floatConstruct(hash(floatBitsToUint(x))); }
float random( vec2  v ) { return floatConstruct(hash(floatBitsToUint(v))); }
float random( vec3  v ) { return floatConstruct(hash(floatBitsToUint(v))); }
float random( vec4  v ) { return floatConstruct(hash(floatBitsToUint(v))); }

// Distance to boundary of screen, working from the inside
float sdf(in vec2 uv)
{
    return min(min(uv.x, uv.y), min(1. - uv.x, 1. - uv.y));
}

// Compute colour at boundary...
vec3 boundary(in vec2 xy)
{
    return xy.x > 0.5 ? vec3(0.5 * xy.y, 0.5 + 0.5 * cos(15.0 * xy.y), 0.0)
                      : vec3(1.0 - xy.y, 0.5 + 0.5 * sin(22.0 * xy.y), 1.0);
}

// ...and interpolate to rest of image.
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec3 total = vec3(0.0, 0.0, 0.0);

    // number_of_walks is essentially the number of samples
    // per pixel so noise is proportional
    // to 1/sqrt(this).
    const int number_of_walks = 200;
    
    const int steps_per_walk = 5; // Seems low but looks fine
    
    // Average n random walks
    for (int walk = 0; walk < number_of_walks; ++walk)
    {
        vec2 uv = fragCoord / iResolution.xy;
        
        // Take m steps...
        for (int step = 0; step < steps_per_walk; ++step)
        {
            // ...each of which has a random direction but whose
            // length is the radius of the largest circle you can
            // place at the current point just touching the
            // boundary. That's more efficient of taking lots
            // of tiny steps and it's obvious that taking a random
            // walk from the centre of a circle hits the boundary
            // at uniformly distributed points.
            float r = sdf(uv);
            float theta = 2. * 3.14159265 * random(uv + vec2(step, walk));
            uv.x += r * cos(theta);
            uv.y += r * sin(theta);
        }
     
        total += boundary(uv);
    }
    
    total = total / float(number_of_walks);

    // Output to screen
    fragColor = vec4(total, 1.0);
}