Skip to content

Instantly share code, notes, and snippets.

@dpiponi
Created May 18, 2022 03:39
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save dpiponi/cb585cfc36462ff679e93c5024717c5d to your computer and use it in GitHub Desktop.
Walk-on-spheres
// 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);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment