Shader Dev

Shader Craft

A unified skill covering 36 GLSL shader techniques (ShaderToy-compatible) for real-time visual effects.

Invocation

/shader-dev <request>

$ARGUMENTS contains the user's request (e.g. "create a raymarched SDF scene with soft shadows").

Skill Structure

shader-dev/
├── SKILL.md                      # Core skill (this file)
├── techniques/                   # Implementation guides (read per routing table)
│   ├── ray-marching.md           # Sphere tracing with SDF
│   ├── sdf-3d.md                 # 3D signed distance functions
│   ├── lighting-model.md         # PBR, Phong, toon shading
│   ├── procedural-noise.md       # Perlin, Simplex, FBM
│   └── ...                       # 34 more technique files
└── reference/                    # Detailed guides (read as needed)
    ├── ray-marching.md           # Math derivations & advanced patterns
    ├── sdf-3d.md                 # Extended SDF theory
    ├── lighting-model.md         # Lighting math deep-dive
    ├── procedural-noise.md       # Noise function theory
    └── ...                       # 34 more reference files

How to Use

1. Read the Technique Routing Table below to identify which technique(s) match the user's request 2. Read the relevant file(s) from techniques/ — each file contains core principles, implementation steps, and complete code templates 3. If you need deeper understanding (math derivations, advanced patterns), follow the reference link at the bottom of each technique file to reference/ 4. Apply the WebGL2 Adaptation Rules below when generating standalone HTML pages

Technique Routing Table

Technique Index

Geometry & SDF

• sdf-2d — 2D signed distance functions for shapes, UI, anti-aliased rendering
• sdf-3d — 3D signed distance functions for real-time implicit surface modeling
• csg-boolean-operations — Constructive solid geometry: union, subtraction, intersection with smooth blending
• domain-repetition — Infinite space repetition, folding, and limited tiling
• domain-warping — Distort domains with noise for organic, flowing shapes
• sdf-tricks — SDF optimization, bounding volumes, binary search refinement, hollowing, layered edges, debug visualization

Ray Casting & Lighting

• ray-marching — Sphere tracing with SDF for 3D scene rendering
• analytic-ray-tracing — Closed-form ray-primitive intersections (sphere, plane, box, torus)
• path-tracing-gi — Monte Carlo path tracing for photorealistic global illumination
• lighting-model — Phong, Blinn-Phong, PBR (Cook-Torrance), and toon shading
• shadow-techniques — Hard shadows, soft shadows (penumbra estimation), cascade shadows
• ambient-occlusion — SDF-based AO, screen-space AO approximation
• normal-estimation — Finite-difference normals, tetrahedron technique

Simulation & Physics

• fluid-simulation — Navier-Stokes fluid solver with advection, diffusion, pressure projection
• simulation-physics — GPU-based physics: springs, cloth, N-body gravity, collision
• particle-system — Stateless and stateful particle systems (fire, rain, sparks, galaxies)
• cellular-automata — Game of Life, reaction-diffusion (Turing patterns), sand simulation

Natural Phenomena

• water-ocean — Gerstner waves, FFT ocean, caustics, underwater fog
• terrain-rendering — Heightfield ray marching, FBM terrain, erosion
• atmospheric-scattering — Rayleigh/Mie scattering, god rays, SSS approximation
• volumetric-rendering — Volume ray marching for clouds, fog, fire, explosions

Procedural Generation

• procedural-noise — Value noise, Perlin, Simplex, Worley, FBM, ridged noise
• procedural-2d-pattern — Brick, hexagon, truchet, Islamic geometric patterns
• voronoi-cellular-noise — Voronoi diagrams, Worley noise, cracked earth, crystal
• fractal-rendering — Mandelbrot, Julia sets, 3D fractals (Mandelbox, Mandelbulb)
• color-palette — Cosine palettes, HSL/HSV/Oklab, dynamic color mapping

Post-Processing & Infrastructure

• post-processing — Bloom, tone mapping (ACES, Reinhard), vignette, chromatic aberration, glitch
• multipass-buffer — Ping-pong FBO setup, state persistence across frames
• texture-sampling — Bilinear, bicubic, mipmap, procedural texture lookup
• matrix-transform — Camera look-at, projection, rotation, orbit controls
• polar-uv-manipulation — Polar/log-polar coordinates, kaleidoscope, spiral mapping
• anti-aliasing — SSAA, SDF analytical AA, temporal anti-aliasing (TAA), FXAA post-process
• camera-effects — Depth of field (thin lens), motion blur, lens distortion, film grain, vignette
• texture-mapping-advanced — Biplanar mapping, texture repetition avoidance, ray differential filtering

Audio

• sound-synthesis — Procedural audio in GLSL: oscillators, envelopes, filters, FM synthesis

Debugging & Validation

• webgl-pitfalls — Common WebGL2/GLSL errors: fragCoord, main() wrapper, function order, macro limitations, uniform null

WebGL2 Adaptation Rules

All technique files use ShaderToy GLSL style. When generating standalone HTML pages, apply these adaptations:

Shader Version & Output

• Use canvas.getContext("webgl2")
• Shader first line: #version 300 es, fragment shader adds precision highp float;
• Fragment shader must declare: out vec4 fragColor;
• Vertex shader: attribute → in, varying → out
• Fragment shader: varying → in, gl_FragColor → fragColor, texture2D() → texture()

Fragment Coordinate

• Use `gl_FragCoord.xy` instead of fragCoord (WebGL2 does not have fragCoord built-in)

// WRONG
vec2 uv = (2.0 * fragCoord - iResolution.xy) / iResolution.y;
// CORRECT
vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / iResolution.y;

main() Wrapper for ShaderToy Templates

• ShaderToy uses void mainImage(out vec4 fragColor, in vec2 fragCoord)
• WebGL2 requires standard void main() entry point — always wrap mainImage:

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    // shader code...
    fragColor = vec4(col, 1.0);
}

void main() {
    mainImage(fragColor, gl_FragCoord.xy);
}

Function Declaration Order

• GLSL requires functions to be declared before use — either declare before use or reorder:

// WRONG — getAtmosphere() calls getSunDirection() before it's defined
vec3 getAtmosphere(vec3 dir) { return getSunDirection(); } // Error!
vec3 getSunDirection() { return normalize(vec3(1.0)); }

// CORRECT — define callee first
vec3 getSunDirection() { return normalize(vec3(1.0)); }
vec3 getAtmosphere(vec3 dir) { return getSunDirection(); } // Works

Macro Limitations

• #define cannot use function calls — use const instead:

// WRONG
#define SUN_DIR normalize(vec3(0.8, 0.4, -0.6))

// CORRECT
const vec3 SUN_DIR = vec3(0.756, 0.378, -0.567); // Pre-computed normalized value

Script Tag Extraction

• When extracting shader source from <script> tags, ensure #version is the first character — use .trim():

const fs = document.getElementById('fs').text.trim();

Common Pitfalls

• Unused uniforms: Compiler may optimize away unused uniforms, causing gl.getUniformLocation() to return null — always use uniforms in a way the compiler cannot optimize out
• Loop indices: Use runtime constants in loops, not #define macros in some ES versions
• Terrain functions: Functions like terrainM(vec2) need XZ components — use terrainM(pos.xz + offset) not terrainM(pos + offset)

HTML Page Setup

When generating a standalone HTML page:

• Canvas fills the entire viewport, auto-resizes on window resize
• Page background black, no scrollbars: body { margin: 0; overflow: hidden; background: #000; }
• Implement ShaderToy-compatible uniforms: iTime, iResolution, iMouse, iFrame
• For multi-pass effects (Buffer A/B), use WebGL2 framebuffer + ping-pong (see multipass-buffer technique)

Common Pitfalls

JS Variable Declaration Order (TDZ — causes white screen crash)

let/const variables must be declared at the top of the <script> block, before any function that references them:

// 1. State variables FIRST
let frameCount = 0;
let startTime = Date.now();

// 2. Canvas/GL init, shader compile, FBO creation
const canvas = document.getElementById('canvas');
const gl = canvas.getContext('webgl2');
// ...

// 3. Functions and event bindings LAST
function resize() { /* can now safely reference frameCount */ }
function render() { /* ... */ }
window.addEventListener('resize', resize);

Reason: let/const have a Temporal Dead Zone — referencing them before declaration throws ReferenceError, causing a white screen.

GLSL Compilation Errors (self-check after writing shaders)

• Function signature mismatch: Call must exactly match definition in parameter count and types. If defined as float fbm(vec3 p), cannot call fbm(uv) with a vec2
• Reserved words as variable names: Do not use: patch, cast, sample, filter, input, output, common, partition, active
• Strict type matching: vec3 x = 1.0 is illegal — use vec3 x = vec3(1.0); cannot use .z to access a vec2
• No ternary on structs: ESSL does not allow ternary operator on struct types — use if/else instead

Performance Budget

Deployment environments may use headless software rendering with limited GPU power. Stay within these limits:

• Ray marching main loop: ≤ 128 steps
• Volume sampling / lighting inner loops: ≤ 32 steps
• FBM octaves: ≤ 6 layers
• Total nested loop iterations per pixel: ≤ 1000 (exceeding this freezes the browser)

Quick Recipes

Common effect combinations — complete rendering pipelines assembled from technique modules.

Photorealistic SDF Scene

1. Geometry: sdf-3d (extended primitives) + csg-boolean-operations (cubic/quartic smin) 2. Rendering: ray-marching + normal-estimation (tetrahedron method) 3. Lighting: lighting-model (outdoor three-light model) + shadow-techniques (improved soft shadow) + ambient-occlusion 4. Atmosphere: atmospheric-scattering (height-based fog with sun tint) 5. Post: post-processing (ACES tone mapping) + anti-aliasing (2x SSAA) + camera-effects (vignette)

Organic / Biological Forms

1. Geometry: sdf-3d (extended primitives + deformation operators: twist, bend) + csg-boolean (gradient-aware smin for material blending) 2. Detail: procedural-noise (FBM with derivatives) + domain-warping 3. Surface: lighting-model (subsurface scattering approximation via half-Lambert)

Procedural Landscape

1. Terrain: terrain-rendering + procedural-noise (erosion FBM with derivatives) 2. Texturing: texture-mapping-advanced (biplanar mapping + no-tile) 3. Sky: atmospheric-scattering (Rayleigh/Mie + height fog) 4. Water: water-ocean (Gerstner waves) + lighting-model (Fresnel reflections)

Stylized 2D Art

1. Shapes: sdf-2d (extended library) + sdf-tricks (layered edges, hollowing) 2. Color: color-palette (cosine palettes) + polar-uv-manipulation (kaleidoscope) 3. Polish: anti-aliasing (SDF analytical AA) + post-processing (bloom, chromatic aberration)

Shader Debugging Techniques

Visual debugging methods — temporarily replace your output to diagnose issues.

SDF Ambient Occlusion — Detailed Reference

This document is a detailed supplement to SKILL.md, containing a complete step-by-step tutorial, mathematical derivations, variant analysis, and advanced usage.

Prerequisites

• GLSL basic syntax (uniform, varying, function definitions)
• Signed Distance Field (SDF) concept: map(p) returns the distance from point p to the nearest surface
• Raymarching basic loop: marching along a ray to find surface intersections
• Surface normal computation: Obtaining the normal direction via SDF gradient (finite differences)
• Vector math fundamentals: dot product, normalization, vector addition/subtraction

Core Principles in Detail

The core idea of SDF ambient occlusion: Sample the SDF at multiple distances along the surface normal and compare the "expected distance" with the "actual distance" to estimate the degree of occlusion.

For a point P on the surface with normal N, at distance h:

• Expected distance = h (if the surroundings are completely open, the SDF value should equal the distance to the surface)
• Actual distance = map(P + N × h) (real SDF value)
• Occlusion contribution = h - map(P + N × h) (the larger the difference, the more nearby geometry is occluding)

The final result is a weighted sum of occlusion contributions from multiple sample points, yielding a [0, 1] occlusion factor:

• 1.0 = no occlusion (bright)
• 0.0 = fully occluded (dark corner)

Key mathematical formula (additive accumulation form):

AO = 1 - k × Σ(weight_i × max(0, h_i - map(P + N × h_i)))

Where weight_i typically decays exponentially or geometrically (closer samples have higher weight), and k is a global intensity coefficient.

Implementation Steps in Detail

Step 1: Build the Base SDF Scene

What: Define a map() function that returns the signed distance value for any point in space.

Why: AO computation relies entirely on SDF queries, so a working distance field is needed first.

float map(vec3 p) {
    float d = p.y; // Ground plane
    d = min(d, length(p - vec3(0.0, 1.0, 0.0)) - 1.0); // Sphere
    d = min(d, length(vec2(length(p.xz) - 1.5, p.y - 0.5)) - 0.4); // Torus
    return d;
}

Step 2: Compute Surface Normal

What: Compute the normal direction via finite difference approximation of the SDF gradient.

Why: AO sampling probes outward along the normal direction; the normal determines the sampling direction.

vec3 calcNormal(vec3 p) {
    vec2 e = vec2(0.001, 0.0);
    return normalize(vec3(
        map(p + e.xyy) - map(p - e.xyy),
        map(p + e.yxy) - map(p - e.yxy),
        map(p + e.yyx) - map(p - e.yyx)
    ));
}

Step 3: Implement Classic Normal-Direction AO (Additive Accumulation)

What: Sample the SDF at 5 distances along the normal direction, accumulating occlusion.

Why: This is a classic method — the most concise and efficient SDF-AO implementation. 5 samples strike an excellent balance between quality and performance. The weight decays at 0.95 exponentially, giving closer samples more influence (near-surface occlusion is more perceptually important).

// Classic AO
float calcAO(vec3 pos, vec3 nor) {
    float occ = 0.0;
    float sca = 1.0; // Initial weight
    for (int i = 0; i < 5; i++) {
        float h = 0.01 + 0.12 * float(i) / 4.0; // Sample distance: 0.01 ~ 0.13
        float d = map(pos + h * nor);             // Actual SDF distance
        occ += (h - d) * sca;                     // Accumulate (expected - actual) × weight
        sca *= 0.95;                              // Weight decay
    }
    return clamp(1.0 - 3.0 * occ, 0.0, 1.0);
}

Step 4: Apply AO to Lighting

What: Multiply the AO factor into ambient and indirect light components.

Why: AO simulates the degree to which indirect light is occluded. Physically, it should only affect ambient/indirect light, not the direct light source's diffuse and specular (direct light occlusion is handled by shadows). However, in practice AO is often multiplied into all lighting for a stronger visual effect.

float ao = calcAO(pos, nor);

// Method A: Affect only ambient light (physically correct)
vec3 ambient = vec3(0.2, 0.3, 0.5) * ao;
vec3 color = diffuse * shadow + ambient;

// Method B: Affect all lighting (stronger visual effect)
vec3 color = (diffuse * shadow + ambient) * ao;

// Method C: Combined with sky visibility bias
float skyVis = 0.5 + 0.5 * nor.y; // Upward-facing surfaces are brighter
vec3 color = diffuse * shadow + ambient * ao * skyVis;

Step 5: Raymarching Main Loop Integration

What: Integrate AO into the complete raymarching pipeline.

Why: AO is part of the lighting computation and needs to be calculated after hitting a surface but before final output.

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    // ... camera setup, ray generation ...

    // Raymarching loop
    float t = 0.0;
    for (int i = 0; i < 128; i++) {
        vec3 p = ro + rd * t;
        float d = map(p);
        if (d < 0.001) break;
        t += d;
        if (t > 100.0) break;
    }

    // Compute lighting on hit
    vec3 col = vec3(0.0);
    if (t < 100.0) {
        vec3 pos = ro + rd * t;
        vec3 nor = calcNormal(pos);
        float ao = calcAO(pos, nor);

        // Lighting
        vec3 lig = normalize(vec3(1.0, 0.8, -0.6));
        float dif = clamp(dot(nor, lig), 0.0, 1.0);
        float sky = 0.5 + 0.5 * nor.y;
        col = vec3(1.0) * dif + vec3(0.2, 0.3, 0.5) * sky * ao;
    }

    fragColor = vec4(col, 1.0);
}

Variant Details

Variant 1: Multiplicative AO

Difference from base version: Starts at 1.0 and progressively multiplies down, rather than using additive accumulation then inverting. The multiplicative form naturally guarantees the result stays in [0, 1], avoids the need for clamping, and provides more natural falloff for multiple overlapping occlusions.

Source: Multiplicative accumulation approach

// Multiplicative AO
float calcAO_multiplicative(vec3 pos, vec3 nor) {
    float ao = 1.0;
    float dist = 0.0;
    for (int i = 0; i <= 5; i++) {
        dist += 0.1; // Uniform step of 0.1
        float d = map(pos + nor * dist);
        ao *= 1.0 - max(0.0, (dist - d) * 0.2 / dist);
    }
    return ao;
}

Variant 2: Multi-Scale AO

Difference from base version: Exponentially increases sampling distances (0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4), computing short-range and long-range occlusion separately. Short-range AO reveals contact shadows and surface detail; long-range AO reveals large-scale environmental occlusion. Fully unrolled with no loops, making it GPU-efficient.

Source: Multi-scale sampling approach

// Multi-scale AO
float calcAO_multiscale(vec3 pos, vec3 nor) {
    // Short-range AO (contact shadows)
    float aoS = 1.0;
    aoS *= clamp(map(pos + nor * 0.1) * 10.0, 0.0, 1.0);  // Adjustable: distance 0.1, weight 10.0
    aoS *= clamp(map(pos + nor * 0.2) * 5.0,  0.0, 1.0);  // Adjustable: distance 0.2, weight 5.0
    aoS *= clamp(map(pos + nor * 0.4) * 2.5,  0.0, 1.0);  // Adjustable: distance 0.4, weight 2.5
    aoS *= clamp(map(pos + nor * 0.8) * 1.25, 0.0, 1.0);  // Adjustable: distance 0.8, weight 1.25

    // Long-range AO (large-scale occlusion)
    float ao = aoS;
    ao *= clamp(map(pos + nor * 1.6) * 0.625,  0.0, 1.0); // Adjustable: distance 1.6
    ao *= clamp(map(pos + nor * 3.2) * 0.3125, 0.0, 1.0); // Adjustable: distance 3.2
    ao *= clamp(map(pos + nor * 6.4) * 0.15625,0.0, 1.0);  // Adjustable: distance 6.4

    return max(0.035, pow(ao, 0.3)); // pow compresses dynamic range, min prevents total black
}

Variant 3: Jittered Sampling AO

Difference from base version: Adds hash-based jitter on top of uniform sample positions, breaking the banding artifacts caused by fixed sample spacing. Also uses a 1/(1+l) distance-decay weight so farther samples have less influence.

Source: Jittered sampling approach

// Jittered sampling AO
float hash(float n) { return fract(sin(n) * 43758.5453); }

float calcAO_jittered(vec3 pos, vec3 nor, float maxDist) {
    float ao = 0.0;
    const float nbIte = 6.0;          // Adjustable: number of samples
    for (float i = 1.0; i < nbIte + 0.5; i++) {
        float l = (i + hash(i)) * 0.5 / nbIte * maxDist; // Jittered sample position
        ao += (l - map(pos + nor * l)) / (1.0 + l);       // Distance-decay weight
    }
    return clamp(1.0 - ao / nbIte, 0.0, 1.0);
}
// Usage example: calcAO_jittered(pos, nor, 4.0)

Variant 4: Hemispherical Random Direction AO

Difference from base version: Instead of sampling only along the normal direction, generates multiple random directions within the normal hemisphere. Closer to the true physical model of ambient occlusion (light arriving from all directions in the hemisphere), but requires more samples (typically 32) for smooth results.

Source: Hemispherical random direction approach

// Hemispherical random direction AO
vec2 hash2(float n) {
    return fract(sin(vec2(n, n + 1.0)) * vec2(43758.5453, 22578.1459));
}

float calcAO_hemisphere(vec3 pos, vec3 nor, float seed) {
    float occ = 0.0;
    for (int i = 0; i < 32; i++) {                              // Adjustable: sample count (16~64)
        float h = 0.01 + 4.0 * pow(float(i) / 31.0, 2.0);      // Quadratic distribution biased toward near-field
        vec2 an = hash2(seed + float(i) * 13.1) * vec2(3.14159, 6.2831); // Random spherical coordinates
        vec3 dir = vec3(sin(an.x) * sin(an.y), sin(an.x) * cos(an.y), cos(an.x));
        dir *= sign(dot(dir, nor));                               // Flip to normal hemisphere
        occ += clamp(5.0 * map(pos + h * dir) / h, -1.0, 1.0); // Signed occlusion
    }
    return clamp(occ / 32.0, 0.0, 1.0);
}

Variant 5: Fibonacci Sphere Uniform Hemisphere AO

Difference from base version: Uses Fibonacci sphere points instead of random directions, achieving quasi-uniform hemisphere sampling distribution. Avoids the clustering problem of pure random sampling, yielding higher quality at the same sample count. Can also be paired with a separate directional occlusion function (e.g., SSS/soft shadow) for multi-level occlusion.

Source: Fibonacci sphere sampling approach

// Fibonacci sphere sampling AO
vec3 forwardSF(float i, float n) {
    const float PI  = 3.141592653589793;
    const float PHI = 1.618033988749895;
    float phi = 2.0 * PI * fract(i / PHI);
    float zi = 1.0 - (2.0 * i + 1.0) / n;
    float sinTheta = sqrt(1.0 - zi * zi);
    return vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, zi);
}

float hash1(float n) { return fract(sin(n) * 43758.5453); }

float calcAO_fibonacci(vec3 pos, vec3 nor) {
    float ao = 0.0;
    for (int i = 0; i < 32; i++) {                         // Adjustable: sample count
        vec3 ap = forwardSF(float(i), 32.0);
        float h = hash1(float(i));
        ap *= sign(dot(ap, nor)) * h * 0.1;                // Flip to hemisphere + random scale
        ao += clamp(map(pos + nor * 0.01 + ap) * 3.0, 0.0, 1.0);
    }
    ao /= 32.0;
    return clamp(ao * 6.0, 0.0, 1.0);
}

Performance Optimization Details

Bottleneck Analysis

The performance bottleneck of SDF-AO lies almost entirely in SDF sample count — each map() call is a full scene distance computation. For complex scenes, this can be very expensive.

Optimization Techniques

1. Reduce Sample Count

Classic normal-direction AO only needs 3~5 samples for acceptable quality. Hemispherical sampling is more physically correct but requires 16~32 samples; use it when the performance budget allows.

2. Early Exit Optimization

Exit the loop early when accumulated occlusion is already large enough, avoiding unnecessary SDF computations.

if (occ > 0.35) break; // Early exit when heavily occluded

3. Unroll Loops

For fixed sample counts (especially 4~7), manually unrolling loops avoids branch overhead and is GPU-friendly. The multi-scale AO variant fully unrolls 7 samples.

4. Simplify AO for Distant Objects

Objects far from the camera can use fewer AO samples or skip AO entirely.

float aoSteps = mix(5.0, 2.0, clamp(t / 50.0, 0.0, 1.0));

5. Precompilation Switches

Use #ifdef to disable AO in debug or low-performance modes.

#ifdef ENABLE_AMBIENT_OCCLUSION
    float ao = calcAO(pos, nor);
#else
    float ao = 1.0;
#endif

6. Hand-Painted Pseudo-AO Blending

For static or semi-static scenes, pseudo-AO values (based on material ID or position) can be precomputed and blended with real-time AO to reduce runtime computation.

float focc = /* preset occlusion based on material */;
float finalAO = calcAO(pos, nor) * focc;

7. SDF Simplification

A simplified version of map() (ignoring small details) can be used for AO sampling, since AO is inherently low-frequency information.

Combination Suggestions in Detail

1. AO + Soft Shadow

The most common combination. AO handles indirect light occlusion (corners, crevices); soft shadows handle direct light occlusion. Simply multiply the two:

float sha = calcShadow(pos, lightDir, 0.02, 20.0, 8.0);
float ao = calcAO(pos, nor);
col = diffuse * sha + ambient * ao; // Each handles its own domain
// Or more simply:
col = lighting * sha * ao;

2. AO + Sky Visibility

Use the normal's y component to estimate the degree of upward-facing, multiplied with AO to simulate sky light occlusion:

float skyVis = 0.5 + 0.5 * nor.y;
col += skyColor * ao * skyVis;

3. AO + Subsurface Scattering / Bounce Light

AO can modulate bounce light and SSS intensity (occluded areas also don't receive bounce light):

float bou = clamp(-nor.y, 0.0, 1.0); // Downward-facing surfaces receive ground bounce
col += bounceColor * bou * ao;
col += sssColor * sss * (0.05 + 0.95 * ao); // SSS also modulated by AO

4. AO + Convexity / Corner Detection

The same SDF probing loop can sample both outward (+N) and inward (-N), yielding AO and convexity information respectively, useful for edge highlights or wear effects:

vec2 aoAndCorner = getOcclusion(pos, nor); // .x = AO, .y = convexity
col *= aoAndCorner.x;                       // AO darkening
col = mix(col, edgeColor, aoAndCorner.y);   // Convexity coloring

5. AO + Fresnel Environment Reflection

AO should also modulate the environment reflection term; otherwise concave areas will show unnatural bright environment reflections:

float fre = pow(1.0 - max(dot(rd, nor), 0.0), 5.0);
col += envColor * fre * ao; // Reduce environment reflection in occluded areas

Analytic Ray Tracing - Detailed Reference

This document is a detailed supplement to SKILL.md, covering prerequisite knowledge, step-by-step tutorial, mathematical derivations, and advanced usage.

Prerequisites

• Vector math fundamentals: Dot product dot(), cross product cross(), vector normalization normalize()
• Quadratic equation solving: Discriminant b²-4ac, meaning of the two roots
• Ray parametric representation: P(t) = ro + t * rd, where ro is the ray origin, rd is the direction, t is the distance
• GLSL fundamentals: struct, inout parameters, vec3/vec4 operations
• ShaderToy framework: mainImage() function, iResolution, iTime, and other uniforms

Use Cases (Complete List)

• When rendering scenes composed of geometric primitives (spheres, planes, boxes, cylinders, tori, etc.)
• When precise surface intersection points, normals, and distances are needed (no iterative approximation required)
• When efficient ray intersection is needed in real-time rendering (several times faster than ray marching)
• Building the underlying geometric engine for ray tracers and path tracers
• Creating visualization effects for hard-surface modeling (jewelry, mechanical parts, chess scenes, etc.)
• Scenes requiring precise shadows, reflections, and refractions (analytic solutions have no sampling error)

Core Principles in Detail

The core idea of analytic ray tracing is: substitute the ray equation P(t) = O + tD into the implicit equation of the geometric body, obtaining an algebraic equation in t, then solve it using closed-form formulas.

Unified Framework

All analytic intersection functions follow the same pattern:

1. Set up equation: Substitute the ray parametric form into the geometry's implicit equation 2. Simplify and solve: Use algebraic identities to reduce to a standard form (quadratic/quartic equation) 3. Discriminant check: Discriminant < 0 indicates no intersection 4. Select nearest intersection: Take the smallest positive root satisfying distance constraints 5. Compute normal: Evaluate the gradient of the implicit equation at the intersection point

Key Mathematical Formulas

Sphere |P-C|² = r² → quadratic equation: t² + 2bt + c = 0

Plane N·P + d = 0 → linear equation: t = -(N·O + d) / (N·D)

Box Intersection of three pairs of parallel planes → Slab Method: tN = max(t1.x, t1.y, t1.z), tF = min(t2.x, t2.y, t2.z)

Ellipsoid |P/R|² = 1 → sphere intersection in scaled space

Torus (|P_xy| - R)² + P_z² = r² → quartic equation, solved via resolvent cubic

Implementation Steps in Detail

Step 1: Ray Generation

What: Generate a ray from the camera position through each pixel.

Why: This is the starting point of ray tracing. Each pixel corresponds to a ray from the camera through the near plane. The standard approach is to construct a camera coordinate system (right, up, forward) and map normalized screen coordinates to world-space directions.

// Construct camera ray
vec3 generateRay(vec2 fragCoord, vec2 resolution, vec3 ro, vec3 ta) {
    vec2 p = (2.0 * fragCoord - resolution) / resolution.y;

    // Build camera coordinate system
    vec3 cw = normalize(ta - ro);               // forward
    vec3 cu = normalize(cross(cw, vec3(0, 1, 0))); // right
    vec3 cv = cross(cu, cw);                    // up

    float fov = 1.5; // Adjustable: field of view control (larger = narrower angle)
    vec3 rd = normalize(p.x * cu + p.y * cv + fov * cw);
    return rd;
}

Step 2: Ray-Sphere Intersection

What: Compute the exact intersection of a ray with a sphere. This is the most fundamental and commonly used intersection function.

Why: Substituting the ray P = O + tD into the sphere equation |P - C|² = r² and expanding yields a quadratic equation in t. The discriminant h = b² - c determines the number of intersections (0, 1, or 2); the smallest positive root is the nearest intersection.

This is a ubiquitous technique, with two common variants:

Code (optimized version, assumes sphere centered at origin):

// Ray-sphere intersection (optimized version for sphere at origin)
// ro: ray origin (sphere center offset already subtracted)
// rd: ray direction (must be normalized)
// r:  sphere radius
// Returns: intersection distance, MAX_DIST if no intersection
float iSphere(vec3 ro, vec3 rd, vec2 distBound, inout vec3 normal, float r) {
    float b = dot(ro, rd);
    float c = dot(ro, ro) - r * r;
    float h = b * b - c;       // Discriminant (optimized: 4a factor omitted)
    if (h < 0.0) return MAX_DIST; // No intersection

    h = sqrt(h);
    float d1 = -b - h;        // Near intersection
    float d2 = -b + h;        // Far intersection

    // Select the nearest intersection within valid range
    if (d1 >= distBound.x && d1 <= distBound.y) {
        normal = normalize(ro + rd * d1);
        return d1;
    } else if (d2 >= distBound.x && d2 <= distBound.y) {
        normal = normalize(ro + rd * d2);
        return d2;
    }
    return MAX_DIST;
}

Code (general version, arbitrary sphere center):

// Ray-sphere intersection (general version, supports arbitrary sphere center)
// sph: vec4(center.xyz, radius)
float sphIntersect(vec3 ro, vec3 rd, vec4 sph) {
    vec3 oc = ro - sph.xyz;
    float b = dot(oc, rd);
    float c = dot(oc, oc) - sph.w * sph.w;
    float h = b * b - c;
    if (h < 0.0) return -1.0;
    return -b - sqrt(h);  // Returns only the near intersection
}

Step 3: Ray-Plane Intersection

What: Compute the intersection of a ray with an infinite plane.

Why: The plane equation N·P + d = 0 substituted with the ray yields a linear equation, solved directly by division. This is the simplest intersection primitive, commonly used for floors, walls, Cornell Boxes, etc. Note: when N·D ≈ 0, the ray is parallel to the plane.

// Ray-plane intersection
// planeNormal: plane normal (must be normalized)
// planeDist:   distance from plane to origin (N·P + planeDist = 0)
float iPlane(vec3 ro, vec3 rd, vec2 distBound, inout vec3 normal,
             vec3 planeNormal, float planeDist) {
    float denom = dot(rd, planeNormal);
    // Only intersects when ray hits the front face of the plane
    if (denom > 0.0) return MAX_DIST;

    float d = -(dot(ro, planeNormal) + planeDist) / denom;

    if (d < distBound.x || d > distBound.y) return MAX_DIST;

    normal = planeNormal;
    return d;
}

// Quick version: horizontal ground plane (y-axis aligned)
float iGroundPlane(vec3 ro, vec3 rd, float height) {
    return -(ro.y - height) / rd.y;
}

Step 4: Ray-Box Intersection (Slab Method)

What: Compute the intersection of a ray with an axis-aligned bounding box (AABB).

Why: The Slab Method treats the box as the intersection of three pairs of parallel planes. It computes the ray's intersection with each pair of planes (tmin, tmax), then takes the maximum of all tmin values and the minimum of all tmax values. If tN > tF or tF < 0, there is no intersection. The normal is determined by which face was hit first.

// Ray-box intersection (Slab Method, optimized version)
// boxSize: box half-size vec3(halfW, halfH, halfD)
float iBox(vec3 ro, vec3 rd, vec2 distBound, inout vec3 normal, vec3 boxSize) {
    vec3 m = sign(rd) / max(abs(rd), 1e-8); // Avoid division by zero
    vec3 n = m * ro;
    vec3 k = abs(m) * boxSize;

    vec3 t1 = -n - k;  // Near plane intersections
    vec3 t2 = -n + k;  // Far plane intersections

    float tN = max(max(t1.x, t1.y), t1.z); // Entry distance into the box
    float tF = min(min(t2.x, t2.y), t2.z); // Exit distance from the box

    if (tN > tF || tF <= 0.0) return MAX_DIST; // No intersection

    if (tN >= distBound.x && tN <= distBound.y) {
        // Normal: determine which face was hit
        normal = -sign(rd) * step(t1.yzx, t1.xyz) * step(t1.zxy, t1.xyz);
        return tN;
    } else if (tF >= distBound.x && tF <= distBound.y) {
        normal = -sign(rd) * step(t1.yzx, t1.xyz) * step(t1.zxy, t1.xyz);
        return tF;
    }
    return MAX_DIST;
}

Step 5: Ray-Ellipsoid Intersection

What: Compute the intersection of a ray with an ellipsoid.

Why: An ellipsoid can be viewed as a sphere scaled differently along each axis. By dividing both the ray origin and direction by the ellipsoid radii R, a sphere intersection is performed in scaled space, then the normal is transformed back to the original space. This "space transformation" technique is one of the core ideas of analytic intersection.

// Ray-ellipsoid intersection
// rad: vec3(rx, ry, rz) three-axis radii
float iEllipsoid(vec3 ro, vec3 rd, vec2 distBound, inout vec3 normal, vec3 rad) {
    // Transform to unit sphere space
    vec3 ocn = ro / rad;
    vec3 rdn = rd / rad;

    float a = dot(rdn, rdn);
    float b = dot(ocn, rdn);
    float c = dot(ocn, ocn);
    float h = b * b - a * (c - 1.0);

    if (h < 0.0) return MAX_DIST;

    float d = (-b - sqrt(h)) / a;

    if (d < distBound.x || d > distBound.y) return MAX_DIST;

    // Normal in original space: gradient of implicit equation |P/R|²=1 → P/(R²)
    normal = normalize((ro + d * rd) / rad);
    return d;
}

Step 6: Ray-Cylinder Intersection

What: Compute the intersection of a ray with a finite cylinder (with end caps).

Why: Cylinder intersection has two parts: (1) project the problem onto a plane perpendicular to the axis, solving a quadratic equation for side surface intersections; (2) check if the intersection is within the finite length, and if not, test the end cap planes.

// Ray-capped cylinder intersection
// pa, pb: two endpoints of the cylinder axis
// ra: cylinder radius
float iCylinder(vec3 ro, vec3 rd, vec2 distBound, inout vec3 normal,
                vec3 pa, vec3 pb, float ra) {
    vec3 ca = pb - pa;          // Cylinder axis vector
    vec3 oc = ro - pa;

    float caca = dot(ca, ca);
    float card = dot(ca, rd);
    float caoc = dot(ca, oc);

    // Project onto plane perpendicular to axis, build quadratic equation
    float a = caca - card * card;
    float b = caca * dot(oc, rd) - caoc * card;
    float c = caca * dot(oc, oc) - caoc * caoc - ra * ra * caca;
    float h = b * b - a * c;

    if (h < 0.0) return MAX_DIST;

    h = sqrt(h);
    float d = (-b - h) / a;

    // Check if side intersection is within finite length
    float y = caoc + d * card;
    if (y > 0.0 && y < caca && d >= distBound.x && d <= distBound.y) {
        normal = (oc + d * rd - ca * y / caca) / ra;
        return d;
    }

    // Test end caps
    d = ((y < 0.0 ? 0.0 : caca) - caoc) / card;
    if (abs(b + a * d) < h && d >= distBound.x && d <= distBound.y) {
        normal = normalize(ca * sign(y) / caca);
        return d;
    }

    return MAX_DIST;
}

Step 7: Scene Intersection & Shading

What: Traverse all objects in the scene, find the nearest intersection, and compute lighting.

Why: Scene traversal in analytic ray tracing is linear — each ray tests all objects sequentially. Through the unified intersection API (distBound parameter), each time a nearer intersection is found, the search range is automatically shortened, achieving implicit culling.

#define MAX_DIST 1e10

// Unified scene intersection function
// Returns vec3(current nearest distance, final intersection distance, material ID)
vec3 worldHit(vec3 ro, vec3 rd, vec2 dist, out vec3 normal) {
    vec3 d = vec3(dist, 0.0); // (distBound.x, distBound.y, matID)
    vec3 tmpNormal;

    // Ground plane
    float t = iPlane(ro, rd, d.xy, normal, vec3(0, 1, 0), 0.0);
    if (t < d.y) { d.y = t; d.z = 1.0; }

    // Sphere
    t = iSphere(ro - vec3(0, 0.5, 0), rd, d.xy, tmpNormal, 0.5);
    if (t < d.y) { d.y = t; d.z = 2.0; normal = tmpNormal; }

    // Box
    t = iBox(ro - vec3(2, 0.5, 0), rd, d.xy, tmpNormal, vec3(0.5));
    if (t < d.y) { d.y = t; d.z = 3.0; normal = tmpNormal; }

    return d;
}

// Basic shading (Lambertian + shadow)
vec3 shade(vec3 pos, vec3 normal, vec3 rd, vec3 albedo) {
    vec3 lightDir = normalize(vec3(-1.0, 0.75, 1.0));

    // Diffuse
    float diff = max(dot(normal, lightDir), 0.0);

    // Ambient
    float amb = 0.5 + 0.5 * normal.y;

    return albedo * (amb * 0.2 + diff * 0.8);
}

Step 8: Reflection & Refraction

What: Implement iterative reflection/refraction for non-recursive ray bounces.

Why: GLSL does not support recursion, so loops are used to simulate multiple bounces. At each bounce, the intersection point plus offset (epsilon) serves as the new ray origin, with the reflected/refracted direction as the new direction. The Fresnel term determines the energy distribution between reflection and refraction.

#define MAX_BOUNCES 4       // Adjustable: number of reflection bounces (more = more realistic but slower)
#define EPSILON 0.001        // Adjustable: self-intersection offset

// Schlick Fresnel approximation
float schlickFresnel(float cosTheta, float F0) {
    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
}

vec3 radiance(vec3 ro, vec3 rd) {
    vec3 color = vec3(0.0);
    vec3 mask = vec3(1.0);
    vec3 normal;

    for (int i = 0; i < MAX_BOUNCES; i++) {
        vec3 res = worldHit(ro, rd, vec2(EPSILON, MAX_DIST), normal);

        if (res.z < 0.5) {
            // No object hit → sky color
            color += mask * vec3(0.6, 0.8, 1.0);
            break;
        }

        vec3 hitPos = ro + rd * res.y;
        vec3 albedo = getAlbedo(res.z);

        // Fresnel reflection coefficient
        float F = schlickFresnel(max(0.0, dot(normal, -rd)), 0.04);

        // Add diffuse contribution
        color += mask * (1.0 - F) * shade(hitPos, normal, rd, albedo);

        // Update mask and ray (reflection)
        mask *= F * albedo;
        rd = reflect(rd, normal);
        ro = hitPos + EPSILON * rd;
    }

    return color;
}

Complete Code Template

For a complete runnable ShaderToy template, see the "Complete Code Template" section in SKILL.md, which includes sphere, plane, and box primitives with support for reflections and Blinn-Phong shading.

The following table describes the adjustable parameters in the template:

Variant Details

Variant 1: Path Tracing

Difference from base version: Replaces deterministic reflection with random hemisphere sampling to achieve global illumination. Requires multi-frame accumulation and random number generation.

Key code:

// Cosine-weighted random hemisphere direction
vec3 cosWeightedRandomHemisphereDirection(vec3 n, inout float seed) {
    vec2 r = hash2(seed);
    vec3 uu = normalize(cross(n, abs(n.y) > 0.5 ? vec3(1,0,0) : vec3(0,1,0)));
    vec3 vv = cross(uu, n);
    float ra = sqrt(r.y);
    float rx = ra * cos(6.2831 * r.x);
    float ry = ra * sin(6.2831 * r.x);
    float rz = sqrt(1.0 - r.y);
    return normalize(rx * uu + ry * vv + rz * n);
}

// Replace reflect in the bounce loop:
rd = cosWeightedRandomHemisphereDirection(normal, seed);
ro = hitPos + EPSILON * rd;
mask *= mat.albedo; // No Fresnel weighting

Variant 2: Analytical Soft Shadow

Difference from base version: Uses the analytical distance from a sphere to the ray to compute soft shadow gradients, without additional sampling.

Key code:

// Sphere soft shadow
float sphSoftShadow(vec3 ro, vec3 rd, vec4 sph) {
    vec3 oc = ro - sph.xyz;
    float b = dot(oc, rd);
    float c = dot(oc, oc) - sph.w * sph.w;
    float h = b * b - c;
    // d: closest distance from ray to sphere surface, t: distance along ray
    float d = sqrt(max(0.0, sph.w * sph.w - h)) - sph.w;
    float t = -b - sqrt(max(h, 0.0));
    return (t > 0.0) ? max(d, 0.0) / t : 1.0;
}

Variant 3: Analytical Antialiasing

Difference from base version: Uses the analytical distance from a sphere to the ray to compute pixel coverage, achieving edge smoothing without multi-sampling.

Key code:

// Sphere distance information (for antialiasing)
vec2 sphDistances(vec3 ro, vec3 rd, vec4 sph) {
    vec3 oc = ro - sph.xyz;
    float b = dot(oc, rd);
    float c = dot(oc, oc) - sph.w * sph.w;
    float h = b * b - c;
    float d = sqrt(max(0.0, sph.w * sph.w - h)) - sph.w; // Closest distance
    return vec2(d, -b - sqrt(max(h, 0.0)));                // (distance, depth)
}

// In rendering, use coverage instead of hard boundary:
float px = 2.0 / iResolution.y; // Pixel size
vec2 dt = sphDistances(ro, rd, sph);
float coverage = 1.0 - clamp(dt.x / (dt.y * px), 0.0, 1.0);
col = mix(bgColor, sphereColor, coverage);

Variant 4: Refraction (with Snell's Law)

Difference from base version: Adds refracted rays; requires detecting whether the ray hits the surface from outside or inside, and flipping the normal accordingly.

Key code:

float refrIndex = 1.5; // Adjustable: index of refraction (glass≈1.5, water≈1.33)

// Add refraction branch in the bounce loop:
bool inside = dot(rd, normal) > 0.0;
vec3 n = inside ? -normal : normal;
float eta = inside ? refrIndex : 1.0 / refrIndex;
vec3 refracted = refract(rd, n, eta);

// Fresnel determines reflection/refraction ratio
float cosI = abs(dot(rd, n));
float F = schlick(cosI, pow((1.0 - eta) / (1.0 + eta), 2.0));

if (refracted != vec3(0.0) && hash1(seed) > F) {
    rd = refracted;
} else {
    rd = reflect(rd, n);
}
ro = hitPos + rd * EPSILON;

Variant 5: Higher-Order Algebraic Surfaces (Quartic Surfaces - Sphere4, Goursat, Torus)

Difference from base version: Substitutes the ray into quartic equations, solving via the resolvent cubic method. Suitable for tori, super-ellipsoids, and similar shapes.

Key code:

// Ray-Sphere4 intersection (|x|⁴+|y|⁴+|z|⁴ = r⁴)
float iSphere4(vec3 ro, vec3 rd, vec2 distBound, inout vec3 normal, float ra) {
    float r2 = ra * ra;
    vec3 d2 = rd*rd, d3 = d2*rd;
    vec3 o2 = ro*ro, o3 = o2*ro;
    float ka = 1.0 / dot(d2, d2);

    float k0 = ka * dot(ro, d3);
    float k1 = ka * dot(o2, d2);
    float k2 = ka * dot(o3, rd);
    float k3 = ka * (dot(o2, o2) - r2 * r2);

    // Reduce to depressed quartic, solve via resolvent cubic
    float c0 = k1 - k0 * k0;
    float c1 = k2 + 2.0 * k0 * (k0 * k0 - 1.5 * k1);
    float c2 = k3 - 3.0 * k0 * (k0 * (k0 * k0 - 2.0 * k1) + 4.0/3.0 * k2);

    float p = c0 * c0 * 3.0 + c2;
    float q = c0 * c0 * c0 - c0 * c2 + c1 * c1;
    float h = q * q - p * p * p * (1.0/27.0);

    if (h < 0.0) return MAX_DIST; // Convex body: only need to handle 2 real roots case

    h = sqrt(h);
    float s = sign(q+h) * pow(abs(q+h), 1.0/3.0);
    float t = sign(q-h) * pow(abs(q-h), 1.0/3.0);

    vec2 v = vec2((s+t) + c0*4.0, (s-t) * sqrt(3.0)) * 0.5;
    float r = length(v);
    float d = -abs(v.y) / sqrt(r + v.x) - c1/r - k0;

    if (d >= distBound.x && d <= distBound.y) {
        vec3 pos = ro + rd * d;
        normal = normalize(pos * pos * pos); // Gradient: 4x³
        return d;
    }
    return MAX_DIST;
}

Performance Optimization Details

1. Distance Bound Pruning

The most important optimization. Each time a nearer intersection is found, distBound.y is shortened, and subsequent objects are automatically skipped:

// distBound.y continuously shrinks with opU
d = opU(d, iSphere(..., d.xy, ...), matId);
d = opU(d, iBox(..., d.xy, ...), matId);   // Automatically skips objects farther than current hit

2. Bounding Sphere / Bounding Box Pre-Test

For complex geometry (tori, Goursat surfaces, etc.), test a simple bounding sphere first to check for possible intersection:

// Test bounding sphere before torus intersection
if (iSphere(ro, rd, distBound, tmpNormal, torus.x + torus.y) > distBound.y) {
    return MAX_DIST; // Bounding sphere missed, skip expensive quartic equation
}

3. Shadow Ray Early Exit

Shadow detection only needs to know "whether there is an occluder," not the nearest intersection, so a simplified intersection function can be used:

// Fast sphere occlusion test (only checks for intersection, no normal computation)
float fastSphIntersect(vec3 ro, vec3 rd, vec3 center, float r) {
    vec3 v = ro - center;
    float b = dot(v, rd);
    float c = dot(v, v) - r * r;
    float d = b * b - c;
    if (d > 0.0) {
        float t = -b - sqrt(d);
        if (t > 0.0) return t;
        t = -b + sqrt(d);
        if (t > 0.0) return t;
    }
    return -1.0;
}

4. Grid Acceleration Structure

For large numbers of identical primitives (e.g., hundreds of spheres), use a spatial grid to accelerate ray traversal:

// 3D DDA grid traversal (for scenes with many spheres)
vec3 pos = floor(ro / GRIDSIZE) * GRIDSIZE;
vec3 ri = 1.0 / rd;
vec3 rs = sign(rd) * GRIDSIZE;
vec3 dis = (pos - ro + 0.5 * GRIDSIZE + rs * 0.5) * ri;

for (int i = 0; i < MAX_STEPS; i++) {
    // Test spheres in current cell
    testSphereInGrid(pos.xz, ro, rd, ...);
    // DDA step to next cell
    vec3 mm = step(dis.xyz, dis.zyx);
    dis += mm * rs * ri;
    pos += mm * rs;
}

5. Avoiding Unnecessary sqrt

Return early when the discriminant is negative, avoiding sqrt() on negative numbers. In some scenarios, the discriminant's sign can be used for coarse pre-filtering:

// Check if ray is heading toward sphere and not inside it
if (c > 0.0 && b > 0.0) return MAX_DIST; // Fast cull

Combination Suggestions in Detail

1. Analytic Intersection + Raymarching SDF

Use analytic primitives for large simple geometry (ground, bounding boxes), and SDF raymarching for complex details (fractals, smooth boolean operations). Analytic intersection provides precise start/end distances, accelerating marching convergence:

float d = iBox(ro, rd, distBound, normal, boxSize); // Analytic box
if (d < MAX_DIST) {
    // Refine with SDF inside the box
    float t = d;
    for (int i = 0; i < 64; i++) {
        float h = sdfScene(ro + t * rd);
        if (h < 0.001) break;
        t += h;
    }
}

2. Analytic Intersection + Volumetric Effects

Use analytic intersection to obtain precise entry/exit distances, then perform volumetric sampling (clouds, fog, subsurface scattering) within that range:

// Use analytic ellipsoid intersection to obtain volume bounds
float tEnter = (-b - sqrt(h)) / a;
float tExit  = (-b + sqrt(h)) / a;
float thickness = tExit - tEnter; // Analytic thickness

// Sample volume within [tEnter, tExit]
vec3 volumeColor = vec3(0.0);
float dt = (tExit - tEnter) / float(VOLUME_STEPS);
for (int i = 0; i < VOLUME_STEPS; i++) {
    vec3 p = ro + rd * (tEnter + float(i) * dt);
    volumeColor += sampleVolume(p) * dt;
}

3. Analytic Intersection + PBR Material System

Analytic intersection provides precise normals and intersection positions, feeding directly into Cook-Torrance and other PBR shading models:

// Cook-Torrance BRDF (requires precise normals)
float D = beckmannDistribution(NdotH, roughness);
float G = geometricAttenuation(NdotV, NdotL, VdotH, NdotH);
float F = fresnelSchlick(VdotH, F0);
vec3 specular = vec3(D * G * F) / (4.0 * NdotV * NdotL);

4. Analytic Intersection + Spatial Transforms

Reuse the same intersection function for transformed geometry by rotating/translating/scaling the ray:

// Rotate object: rotate the ray instead of the object
vec3 localRo = rotateY(ro - objectPos, angle);
vec3 localRd = rotateY(rd, angle);
float t = iBox(localRo, localRd, distBound, localNormal, boxSize);
// Transform normal back to world space
normal = rotateY(localNormal, -angle);

5. Analytic Intersection + Analytical AO / Soft Shadow / Antialiasing

A fully analytic rendering pipeline: intersection, shadows, occlusion, and edge smoothing all use closed-form formulas, producing zero noise:

// Fully analytic pipeline (no random sampling, no noise)
float t = sphIntersect(ro, rd, sph);        // Analytic intersection
float shadow = sphSoftShadow(hitPos, ld, sph); // Analytic soft shadow
float ao = sphOcclusion(hitPos, normal, sph);  // Analytic ambient occlusion
float coverage = sphAntiAlias(ro, rd, sph, px); // Analytic antialiasing

Anti-Aliasing Detailed Reference

Prerequisites

• Understanding of screen-space derivatives (dFdx, dFdy, fwidth)
• Multipass buffer setup (for TAA)
• Basic signal processing concepts

Sampling Theory (Nyquist)

The Nyquist-Shannon theorem states: to accurately represent a signal, sampling rate must be ≥ 2× the highest frequency present. In shader terms:

• Pixel grid = sampling rate
• Procedural detail / edge sharpness = signal frequency
• When detail frequency > pixel frequency → aliasing (moiré, crawling edges)

Solutions: either increase sampling rate (SSAA) or reduce signal frequency (analytical AA, filtering).

SSAA Implementation Details

Jitter Patterns

• Grid: offset = vec2(m, n) / AA - 0.5 — simple, uniform coverage
• Rotated grid (RGSS): 4 samples at rotated positions — better edge coverage for near-horizontal/vertical lines
• Halton sequence: quasi-random low-discrepancy — best coverage for high sample counts

Performance

AA=2 (4 samples) is the practical limit for real-time SDF scenes. AA=3 (9 samples) for offline/screenshot quality only.

SDF Analytical AA Deep Dive

Why fwidth Works

fwidth(d) = abs(dFdx(d)) + abs(dFdy(d)) approximates how much the SDF value changes across one pixel. Using this as the smoothstep width:

• Edge transition spans exactly ~1 pixel regardless of zoom level
• No texture sampling needed — purely analytical
• Works for any SDF shape

Signed Distance to Coverage

For a 2D SDF with value d at a pixel center:

coverage ≈ clamp(0.5 - d / fwidth(d), 0.0, 1.0)

This maps the signed distance to an approximate pixel coverage, equivalent to a box filter over the pixel footprint.

TAA with Neighborhood Clamping

Full TAA pipeline: 1. Jitter: offset pixel center by Halton(2,3) sequence each frame 2. Render: full scene at jittered position → Buffer A 3. Reproject: use motion vectors to find previous frame's pixel for current position 4. Clamp: restrict history color to the min/max of current frame's 3×3 neighborhood (prevents ghosting) 5. Blend: output = mix(current, clampedHistory, 0.9)

Neighborhood Clamping

vec3 minCol = vec3(1e10), maxCol = vec3(-1e10);
for (int x = -1; x <= 1; x++)
for (int y = -1; y <= 1; y++) {
    vec3 s = texelFetch(currentBuffer, ivec2(fragCoord) + ivec2(x,y), 0).rgb;
    minCol = min(minCol, s);
    maxCol = max(maxCol, s);
}
vec3 clampedHistory = clamp(history, minCol, maxCol);

FXAA Algorithm Walkthrough

1. Luma computation: Convert 5 samples (center + NSEW) to luminance 2. Edge detection: lumaRange = lumaMax - lumaMin — skip if below threshold 3. Edge orientation: Compare horizontal vs vertical luma gradients to determine edge direction 4. Sub-pixel blending: Sample along the edge direction at 1/3 and 2/3 offsets 5. Quality: The simplified version uses 2 taps; full FXAA 3.11 uses up to 12 taps along the edge for better endpoint detection

Atmospheric & Subsurface Scattering — Detailed Reference

This document is a detailed supplement to SKILL.md, covering prerequisites, step-by-step explanations, mathematical derivations, variant details, and complete combination code examples.

Prerequisites

Foundational concepts required before using this Skill:

• GLSL Fundamentals: uniforms, varyings, built-in functions
• Vector Math: dot product, cross product, vector normalization
• Ray-Sphere Intersection: given a ray origin and direction, find the intersection distances with a sphere surface
• Physical Meaning of Exponential Functions (Beer-Lambert Law): light attenuates exponentially through a medium, I = I₀ × e^(-σ×d), where σ is the extinction coefficient and d is the distance
• Basic Ray Marching Concepts: advancing step by step along a ray direction, accumulating information at each sample point

Core Principles

Atmospheric scattering simulates the process of photons passing through the atmosphere and colliding with gas molecules/aerosol particles, changing direction. There are three core physical mechanisms:

1. Rayleigh Scattering (Molecular Scattering)

Caused by particles much smaller than the wavelength of light (nitrogen, oxygen molecules). Short wavelengths (blue light) scatter much more strongly than long wavelengths (red light) — this is why the sky is blue and sunsets are red.

The scattering coefficient is inversely proportional to the fourth power of wavelength:

β_R(λ) ∝ 1/λ⁴

Typical sea-level values for Earth: β_R = vec3(5.5e-6, 13.0e-6, 22.4e-6) (RGB channels, in m⁻¹)

Rayleigh Phase Function (describes the angular distribution of light scattering, symmetric front-to-back):

P_R(θ) = 3/(16π) × (1 + cos²θ)

2. Mie Scattering (Aerosol Scattering)

Caused by particles roughly the same size as the wavelength of light (water droplets, dust). Wavelength-independent (all colors scatter equally), but with strong forward scattering characteristics, forming the halo around the sun.

Typical sea-level values for Earth: β_M = vec3(21e-6) (same for all channels)

Henyey-Greenstein Phase Function (describes the strong forward scattering of Mie scattering):

P_HG(θ, g) = (1 - g²) / (4π × (1 + g² - 2g·cosθ)^(3/2))

Where g ∈ (-1, 1) controls forward scattering strength; typical Earth atmosphere value g ≈ 0.76 ~ 0.88.

3. Beer-Lambert Attenuation

Exponential attenuation of light through a medium:

T(A→B) = exp(-∫ σ_e(s) ds)   // Transmittance from A to B

Where σ_e is the extinction coefficient (extinction = scattering + absorption).

Overall Algorithm Flow

March along the view direction (ray march), at each sample point: 1. Compute the atmospheric density at that point (decreases exponentially with altitude) 2. Perform a second march toward the light source to compute the optical depth from the sun to that point 3. Use Beer-Lambert to calculate the sun light intensity reaching that point 4. Use the phase function to compute the amount of light scattered toward the camera 5. Accumulate contributions from all sample points

Implementation Steps

Step 1: Ray-Sphere Intersection

What: Compute the intersection points of the view ray with the atmospheric shell to determine the ray march start/end range.

Why: The atmosphere is a spherical shell around the planet; we only integrate within the shell.

// Ray-sphere intersection, returns distances to two intersection points (t_near, t_far)
// p: ray origin (relative to sphere center), dir: ray direction, r: sphere radius
vec2 raySphereIntersect(vec3 p, vec3 dir, float r) {
    float b = dot(p, dir);
    float c = dot(p, p) - r * r;
    float d = b * b - c;
    if (d < 0.0) return vec2(1e5, -1e5); // No intersection
    d = sqrt(d);
    return vec2(-b - d, -b + d);
}

Derivation: sphere equation |p + t·dir|² = r² expands to t² + 2t·dot(p,dir) + dot(p,p) - r² = 0. Since dir is normalized, a=1 can be omitted, and the two t values are solved directly with the quadratic formula.

Step 2: Define Atmospheric Physical Constants

What: Set the scale parameters and scattering coefficients for the planet and atmosphere.

Why: These physical constants determine the sky's color characteristics. The different RGB values in Rayleigh produce the blue sky (blue channel has the largest scattering coefficient); Mie's uniform values produce white halos (all wavelengths scatter equally).

#define PLANET_RADIUS 6371e3          // Earth radius (m)
#define ATMOS_RADIUS  6471e3          // Atmosphere outer radius (m), about 100km above Earth's radius
#define PLANET_CENTER vec3(0.0)       // Planet center position

// Scattering coefficients (m⁻¹), sea-level values
#define BETA_RAY vec3(5.5e-6, 13.0e-6, 22.4e-6) // Tunable: Rayleigh scattering, changes sky base color
#define BETA_MIE vec3(21e-6)                      // Tunable: Mie scattering, changes halo intensity
#define BETA_OZONE vec3(2.04e-5, 4.97e-5, 1.95e-6) // Tunable: ozone absorption, affects zenith deep blue

// Mie phase function anisotropy parameter
#define MIE_G 0.76   // Tunable: 0.76~0.88, larger = more concentrated sun halo

// Scale heights (m): altitude at which density drops to 1/e
#define H_RAY 8000.0  // Tunable: Rayleigh scale height, larger = thicker atmosphere
#define H_MIE 1200.0  // Tunable: Mie scale height, larger = higher haze layer

// Ozone parameters (optional)
#define H_OZONE 30e3         // Ozone peak altitude
#define OZONE_FALLOFF 4e3    // Ozone falloff width

// Sample step counts
#define PRIMARY_STEPS 32 // Tunable: primary ray steps, more = higher quality
#define LIGHT_STEPS 8    // Tunable: light direction steps

Parameter tuning guide:

• Increase overall BETA_RAY → more vivid sky color
• Modify BETA_RAY RGB ratios → change sky base hue (e.g., increasing the red component produces a more purple sky)
• Increase BETA_MIE → brighter halo around the sun, more haze
• Increase MIE_G → halo more concentrated toward the sun direction (narrower disk)
• Increase H_RAY → effective atmosphere thickness increases, sky color more uniform
• Increase H_MIE → haze layer higher, low-altitude fog effect weakened

Step 3: Implement Phase Functions

What: Compute the probability distribution of light being scattered at different angles.

Why: The Rayleigh phase is symmetrically distributed (scatters both forward and backward); the Mie phase is strongly biased forward. This determines the brightness distribution across the sky — brighter facing the sun (Mie dominant), with some brightness away from the sun (Rayleigh dominant).

// Rayleigh phase function: symmetric front-to-back
float phaseRayleigh(float cosTheta) {
    return 3.0 / (16.0 * 3.14159265) * (1.0 + cosTheta * cosTheta);
}

// Henyey-Greenstein phase function: forward scattering
// g: anisotropy parameter, 0 = isotropic, close to 1 = strong forward scattering
float phaseMie(float cosTheta, float g) {
    float gg = g * g;
    float num = (1.0 - gg) * (1.0 + cosTheta * cosTheta);
    float denom = (2.0 + gg) * pow(1.0 + gg - 2.0 * g * cosTheta, 1.5);
    return 3.0 / (8.0 * 3.14159265) * num / denom;
}

Note: the Mie phase function here uses the Cornette-Shanks improved version (with an additional (1 + cos²θ) term in the numerator and (2 + g²) normalization correction in the denominator), which is more physically accurate than the original HG.

Step 4: Atmospheric Density Sampling

What: Compute the atmospheric particle density at a given point based on altitude.

Why: Atmospheric density decreases exponentially with altitude, and different components (Rayleigh, Mie, ozone) have different decay rates. Rayleigh particles (gas molecules) have a scale height of about 8km, Mie particles (aerosols) are concentrated in the lower layer with a scale height of about 1.2km, and ozone peaks at approximately 30km altitude.

// Returns vec3(rayleigh_density, mie_density, ozone_density)
vec3 atmosphereDensity(vec3 pos, float planetRadius) {
    float height = length(pos) - planetRadius;

    float densityRay = exp(-height / H_RAY);
    float densityMie = exp(-height / H_MIE);

    // Ozone: peaks at ~30km altitude, approximated with Lorentzian distribution
    float denom = (H_OZONE - height) / OZONE_FALLOFF;
    float densityOzone = (1.0 / (denom * denom + 1.0)) * densityRay;

    return vec3(densityRay, densityMie, densityOzone);
}

Mathematical explanation of ozone distribution: 1/(x² + 1) is the form of a Lorentzian/Cauchy distribution, reaching its maximum value of 1 at x=0 (i.e., height = H_OZONE), then symmetrically decaying on both sides. Multiplying by densityRay accounts for ozone also being affected by the overall atmospheric density decrease.

Step 5: Light Direction Optical Depth

What: From a sample point on the primary ray, march toward the sun to the atmosphere edge, accumulating optical depth.

Why: This determines how much the sunlight has been attenuated before reaching that point. At sunset, the light path passes through more atmosphere, and blue light is scattered away (because Rayleigh scattering coefficient's blue component is largest), leaving only red light — this is the physical reason sunsets are red.

// Compute optical depth from pos along sunDir to the atmosphere edge
vec3 lightOpticalDepth(vec3 pos, vec3 sunDir) {
    float atmoDist = raySphereIntersect(pos - PLANET_CENTER, sunDir, ATMOS_RADIUS).y;
    float stepSize = atmoDist / float(LIGHT_STEPS);
    float rayPos = stepSize * 0.5;

    vec3 optDepth = vec3(0.0); // (ray, mie, ozone)

    for (int i = 0; i < LIGHT_STEPS; i++) {
        vec3 samplePos = pos + sunDir * rayPos;
        float height = length(samplePos - PLANET_CENTER) - PLANET_RADIUS;

        // If sample point is below the surface, it's occluded by the planet
        if (height < 0.0) return vec3(1e10); // Fully occluded

        vec3 density = atmosphereDensity(samplePos, PLANET_RADIUS);
        optDepth += density * stepSize;

        rayPos += stepSize;
    }
    return optDepth;
}

stepSize * 0.5 as the starting offset is the midpoint sampling rule, which approximates the integral more accurately than endpoint sampling.

Step 6: Primary Scattering Integral (Core Loop)

What: Ray march along the view direction, computing the in-scattering contribution at each sample point and accumulating.

Why: This is the core of the entire algorithm — integrating all scattered light along the view direction that reaches the eye. Each point's contribution = sunlight reaching that point × density at that point × attenuation from that point to the camera.

Mathematical expression:

L(camera) = ∫[tStart→tEnd] sunIntensity × T(sun→s) × σ_s(s) × P(θ) × T(s→camera) ds

Where T is transmittance, σ_s is the scattering coefficient, and P is the phase function.

vec3 calculateScattering(
    vec3 rayOrigin,    // Camera position
    vec3 rayDir,       // View direction
    float maxDist,     // Maximum distance (scene occlusion)
    vec3 sunDir,       // Sun direction
    vec3 sunIntensity  // Sun intensity
) {
    // Compute ray-atmosphere intersection
    vec2 atmoHit = raySphereIntersect(rayOrigin - PLANET_CENTER, rayDir, ATMOS_RADIUS);
    if (atmoHit.x > atmoHit.y) return vec3(0.0); // Missed atmosphere

    // Compute ray-planet intersection (ground occlusion)
    vec2 planetHit = raySphereIntersect(rayOrigin - PLANET_CENTER, rayDir, PLANET_RADIUS);

    // Determine march range
    float tStart = max(atmoHit.x, 0.0);
    float tEnd = atmoHit.y;
    if (planetHit.x > 0.0) tEnd = min(tEnd, planetHit.x); // Ground occlusion
    tEnd = min(tEnd, maxDist); // Scene object occlusion

    float stepSize = (tEnd - tStart) / float(PRIMARY_STEPS);

    // Precompute phase functions (view-sun angle is constant along the entire ray)
    float cosTheta = dot(rayDir, sunDir);
    float phaseR = phaseRayleigh(cosTheta);
    float phaseM = phaseMie(cosTheta, MIE_G);

    // Accumulators
    vec3 totalRay = vec3(0.0); // Rayleigh in-scatter
    vec3 totalMie = vec3(0.0); // Mie in-scatter
    vec3 optDepthI = vec3(0.0); // View direction optical depth (ray, mie, ozone)

    float rayPos = tStart + stepSize * 0.5;

    for (int i = 0; i < PRIMARY_STEPS; i++) {
        vec3 samplePos = rayOrigin + rayDir * rayPos;

        // 1. Sample density
        vec3 density = atmosphereDensity(samplePos, PLANET_RADIUS) * stepSize;
        optDepthI += density;

        // 2. Compute light direction optical depth
        vec3 optDepthL = lightOpticalDepth(samplePos, sunDir);

        // 3. Beer-Lambert attenuation: total attenuation from sun through this point to camera
        vec3 tau = BETA_RAY * (optDepthI.x + optDepthL.x)
                 + BETA_MIE * 1.1 * (optDepthI.y + optDepthL.y) // 1.1 is Mie extinction/scattering ratio
                 + BETA_OZONE * (optDepthI.z + optDepthL.z);
        vec3 attenuation = exp(-tau);

        // 4. Accumulate in-scattering
        totalRay += density.x * attenuation;
        totalMie += density.y * attenuation;

        rayPos += stepSize;
    }

    // 5. Final color = scattering coefficient × phase function × accumulated scattering
    return sunIntensity * (
        totalRay * BETA_RAY * phaseR +
        totalMie * BETA_MIE * phaseM
    );
}

Key detail explanations:

• 1.1 is the Mie extinction/scattering ratio: Mie particles not only scatter light but also absorb a small amount, so the extinction coefficient ≈ 1.1 × scattering coefficient
• optDepthI records all three components simultaneously for correctly compositing all extinction contributions in the attenuation calculation
• Phase functions are precomputed outside the loop because the angle between view and sun directions is constant along the entire ray

Step 7: Tone Mapping and Output

What: Apply tone mapping and gamma correction to the HDR scattering results.

Why: The scattering calculation outputs HDR linear values (potentially much greater than 1.0), which must be mapped to [0,1] for display. Different tonemapping methods affect the final look:

• Exposure mapping `1 - exp(-x)`: simplest, naturally saturates and never overexposes, but limited highlight detail
• Reinhard: preserves more highlight detail, suitable for high dynamic range scenes
• ACES: cinematic tone mapping, richer colors but more complex implementation

// Method 1: Simple exposure mapping (most common)
vec3 tonemapExposure(vec3 color) {
    return 1.0 - exp(-color); // Natural saturation, never overexposes
}

// Method 2: Reinhard (preserves more highlight detail)
vec3 tonemapReinhard(vec3 color) {
    float l = dot(color, vec3(0.2126, 0.7152, 0.0722));
    vec3 tc = color / (color + 1.0);
    return mix(color / (l + 1.0), tc, tc);
}

// Gamma correction
vec3 gammaCorrect(vec3 color) {
    return pow(color, vec3(1.0 / 2.2));
}

Reinhard implementation detail: uses a blend of luminance l (perceptually weighted) and per-channel mapping tc, balancing color fidelity and highlight detail.

Variant Details

Variant 1: Non-Physical Analytical Approximation (No Ray March)

Difference from the base version: No ray marching at all — uses analytical functions to simulate sky color with extremely high performance. Not based on physical scattering equations, but uses empirical formulas to simulate visual effects.

Use cases: Mobile platforms, backgrounds, scenes with low physical accuracy requirements.

How it works:

• zenithDensity simulates atmospheric density variation with viewing angle (denser looking toward the horizon)
• getSkyAbsorption uses exp2 to simulate atmospheric absorption (similar to Beer-Lambert)
• getMie uses distance falloff + smoothstep to simulate the sun halo
• The final blend considers the sun altitude's effect on the overall sky color tone

Performance comparison: No loops, no ray march — only a small amount of math per pixel, 10-50x faster than the base version.

Variant 2: With Ozone Absorption Layer

Difference from the base version: Adds ozone absorption as a third component, making the zenith deeper blue and introducing subtle purple tones at sunset.

Use cases: Pursuing more physically accurate sky colors.

Physical principle: Ozone primarily absorbs in the Chappuis band (500-700nm, i.e., green and red), which makes the zenith direction (short light path, remaining light after Rayleigh scattering is filtered by ozone) appear deeper blue. At sunset, the long light path makes ozone absorption more significant — after red is Rayleigh-scattered and green is ozone-absorbed, only blue-purple tones remain.

Key modification: Set BETA_OZONE to a non-zero value in the complete template to enable — already built-in.

Variant 3: Subsurface Scattering (SSS)

Difference from the base version: Scatters inside a semi-transparent object rather than in the atmosphere. Estimates object thickness via SDF and controls light transmission with thickness.

Use cases: Candles, skin, jelly, leaves, and other translucent materials.

How it works: 1. Use Snell's law (refract) to calculate the refracted direction after light enters the object 2. March along the refracted direction in the SDF, accumulating negative distance values (SDF is negative inside the object) 3. Greater accumulated negative value means a thicker object, less light transmission 4. Use a power function to control the attenuation curve (pow parameter is tunable)

Tunable parameters:

• IOR (index of refraction): 1.3 (water) ~ 1.5 (glass) ~ 2.0 (gemstone), affects refraction angle
• MAX_SCATTER: maximum scatter march distance, affects SSS penetration depth
• SCATTER_STRENGTH: scattering intensity multiplier
• Step size 0.2: smaller = more accurate but slower

Usage:

float ss = max(0.0, subsurface(hitPos, viewDir, normal));
vec3 sssColor = albedo * smoothstep(0.0, 2.0, pow(ss, 0.6));
finalColor = mix(lambertian, sssColor, 0.7) + specular;

Variant 4: LUT Precomputation Pipeline (Production-Grade)

Difference from the base version: Precomputes Transmittance, Multiple Scattering, and Sky-View into separate LUT textures; at runtime only performs lookups, with extremely high frame rates.

Use cases: Production-grade sky rendering in game engines and real-time applications requiring high frame rates.

Architecture details:

• Buffer A (Transmittance LUT): 256x64 texture, parameterized by (sunCosZenith, height), storing transmittance from a certain height along a direction to the atmosphere edge. This is the most fundamental LUT; all other LUTs depend on it.

• Buffer B (Multiple Scattering LUT): 32x32 texture, precomputing multiple scattering contributions. Single scattering is not accurate enough — in the real atmosphere, light is scattered multiple times. This LUT uses an iterative method to approximate the cumulative effect of multiple scattering.

• Buffer C (Sky-View LUT): 200x200 texture, storing sky colors for all directions. Uses nonlinear height mapping to allocate more precision to the horizon region (where color changes are most dramatic).

• Image Pass: Only looks up the Sky-View LUT + overlays the sun disk; each pixel requires only one texture query.

// Transmittance LUT query (from Hillaire 2020 implementation)
vec3 getValFromTLUT(sampler2D tex, vec2 bufferRes, vec3 pos, vec3 sunDir) {
    float height = length(pos);
    vec3 up = pos / height;
    float sunCosZenithAngle = dot(sunDir, up);
    vec2 uv = vec2(
        256.0 * clamp(0.5 + 0.5 * sunCosZenithAngle, 0.0, 1.0),
        64.0 * max(0.0, min(1.0, (height - groundRadiusMM) / (atmosphereRadiusMM - groundRadiusMM)))
    );
    uv /= bufferRes;
    return texture(tex, uv).rgb;
}

Performance: The Image Pass is nearly O(1); all heavy computation is done in low-resolution LUTs. LUTs can be incrementally updated as the sun angle changes.

Variant 5: Analytical Fast Atmosphere (No Ray March but Supports Aerial Perspective)

Difference from the base version: Uses analytical exponential approximations instead of ray marching, while supporting distance-attenuated aerial perspective effects.

Use cases: Game scenes requiring atmospheric perspective without per-pixel ray marching.

How it works:

• getRayleighMie uses 1 - exp(-x) form to approximate the scattering integral (analytical solution based on Beer-Lambert)
• getLightTransmittance uses multiple exponential term superposition to approximate optical depth at different sun altitudes
• No loops required — only a fixed number of math operations per pixel

// Based on Felix Westin's Fast Atmosphere
void getRayleighMie(float opticalDepth, float densityR, float densityM, out vec3 R, out vec3 M) {
    vec3 C_RAYLEIGH = vec3(5.802, 13.558, 33.100) * 1e-6;
    vec3 C_MIE = vec3(3.996e-6);
    R = (1.0 - exp(-opticalDepth * densityR * C_RAYLEIGH / 2.5)) * 2.5;
    M = (1.0 - exp(-opticalDepth * densityM * C_MIE / 0.5)) * 0.5;
}

// Analytical approximation of light transmittance (replaces ray march)
vec3 getLightTransmittance(vec3 lightDir) {
    vec3 C_RAYLEIGH = vec3(5.802, 13.558, 33.100) * 1e-6;
    vec3 C_MIE = vec3(3.996e-6);
    vec3 C_OZONE = vec3(0.650, 1.881, 0.085) * 1e-6;
    float extinction = exp(-clamp(lightDir.y + 0.05, 0.0, 1.0) * 40.0)
                     + exp(-clamp(lightDir.y + 0.5, 0.0, 1.0) * 5.0) * 0.4
                     + pow(clamp(1.0 - lightDir.y, 0.0, 1.0), 2.0) * 0.02
                     + 0.002;
    return exp(-(C_RAYLEIGH + C_MIE + C_OZONE) * extinction * 1e6);
}

Mathematical basis of the analytical approximation: Treating the atmosphere as a single uniform layer, the scattering integral ∫ e^(-σx) dx has the analytical solution (1 - e^(-σL)) / σ. The 2.5 and 0.5 in the code are empirical scaling factors to make the analytical result visually approximate a full ray march.

Performance Optimization Details

Bottleneck 1: Nested Ray March (O(N×M) Samples)

N primary ray steps × M light direction steps per step = N×M density calculations.

Optimization approaches:

• Reduce step counts: Use PRIMARY_STEPS=12, LIGHT_STEPS=4 on mobile; visual difference is small but performance improvement is significant
• Analytical approximation: Replace the light direction ray march with the Fast Atmosphere approach, reducing complexity from O(N×M) to O(N)
• Transmittance LUT: After precomputation, runtime only performs lookups, reducing complexity to O(N) or even O(1)

Bottleneck 2: Dense exp() and pow() Calls

Multiple exponential function calls at each sample point — these are relatively expensive operations on GPUs.

Optimization approaches:

• Replace Henyey-Greenstein phase function with Schlick approximation:

// Schlick approximation, only 1 division, no pow
float k = 1.55 * g - 0.55 * g * g * g;
float phaseSchlick = (1.0 - k * k) / (4.0 * PI * pow(1.0 + k * cosTheta, 2.0));

• Combine multiple exp calls: exp(a) * exp(b) = exp(a+b), reducing exp call count
• Use exp2 instead of exp in scenarios with lower precision requirements (exp2 is faster on some GPUs)

Bottleneck 3: Full-Screen Per-Pixel Computation

Each pixel independently computes the full scattering.

Optimization approaches:

• Sky-View LUT: Render the sky to a low-resolution LUT (e.g., 200x200), then look up at full resolution. Allocate more resolution near the horizon (nonlinear mapping)
• Half-resolution rendering: Compute scattering at half resolution, then bilinearly upsample. For sky — a low-frequency signal — quality loss is minimal

Bottleneck 4: High Sample Count Needed to Avoid Banding

Low step counts lead to visible banding artifacts.

Optimization approaches:

• Non-uniform stepping: newT = ((i + 0.3) / numSteps) * tMax, offset by 0.3 instead of 0.5 to reduce visual artifacts
• Jittered start offset: startOffset += hash(fragCoord) * stepSize, randomly offsetting the march start per pixel
• Temporal blue noise dithering: Use temporal blue noise to jitter sample positions across frames; combined with TAA, banding is nearly eliminated

Combination Suggestions

1. Atmospheric Scattering + Volumetric Clouds

Atmospheric scattering provides sky background color and light source color; volumetric cloud lighting uses the atmospheric transmittance to determine the sun light color reaching the cloud layer.

Key integration points:

• Setting the maxDist parameter of the atmospheric scattering function to the cloud layer distance achieves correct pre-cloud atmospheric effects
• During cloud layer rendering, use the transmittance LUT to get the sun light color upon reaching the cloud layer
• Sky color behind clouds should be the full atmospheric scattering result

// Pseudo-code example
float cloudDist = rayMarchClouds(rayOrigin, rayDir);
vec3 cloudColor = calculateCloudLighting(cloudPos, sunDir, transmittance);
vec3 skyBehind = calculateScattering(rayOrigin, rayDir, 1e12, sunDir, sunIntensity);
vec3 skyBeforeCloud = calculateScattering(rayOrigin, rayDir, cloudDist, sunDir, sunIntensity);

// Compositing: pre-cloud atmosphere + cloud × cloud opacity + post-cloud sky × transmittance
vec3 final = skyBeforeCloud + cloudColor * cloudAlpha + skyBehind * (1.0 - cloudAlpha) * atmosphereTransmittance;

2. Atmospheric Scattering + SDF Scene

Pass the SDF ray march hit distance as the maxDist parameter to calculateScattering(), and the scene color as sceneColor, to automatically get aerial perspective effects.

// SDF ray march yields hit information
float hitDist = sdfRayMarch(rayOrigin, rayDir);
vec3 sceneColor = shadeSurface(hitPos, normal, lightDir);

// Atmospheric scattering automatically handles perspective
vec3 final = calculateScattering(
    rayOrigin, rayDir, hitDist,
    sceneColor, sunDir, SUN_INTENSITY
);

3. Atmospheric Scattering + God Rays

Adding an occlusion parameter in the scattering integral (via shadow map or additional ray march for occlusion detection) can produce volumetric light beam effects.

// Add occlusion detection in the main loop
for (int i = 0; i < PRIMARY_STEPS; i++) {
    // ... density sampling ...

    // God rays: check if sample point is occluded
    float occlusion = 1.0;
    if (sdfScene(samplePos + sunDir * 0.1) < 0.0) {
        occlusion = 0.0; // Occluded by scene object, no in-scattering
    }

    totalRay += density.x * attenuation * occlusion;
    totalMie += density.y * attenuation * occlusion;
}

The Fast Atmosphere example implements this functionality through the occlusion parameter.

4. Atmospheric Scattering + Terrain Rendering

Use aerial perspective: distant terrain colors blend into atmospheric scattering color based on distance.

Key formula:

// Basic aerial perspective
vec3 finalColor = terrainColor * transmittance + inscattering;

// transmittance: atmospheric transmittance from camera to terrain point
// inscattering: scattered light between camera and terrain point
// Distant objects: transmittance → 0, inscattering dominates → appears blue/gray

5. SSS + PBR Materials

Combine subsurface scattering with GGX microsurface specular and Fresnel reflection. SSS contribution replaces part of the diffuse (via mix), with the specular layer added on top:

// Complete PBR + SSS shading
float fresnel = pow(max(0.0, 1.0 + dot(normal, viewDir)), 5.0);
vec3 diffuse = mix(lambert, sssContribution, 0.7);  // SSS replaces part of diffuse
vec3 final = ambient + albedo * diffuse + specular + fresnel * envColor;

Layering logic: 1. Bottom layer: ambient light 2. Diffuse layer: blend of Lambert and SSS (SSS allows light to pass through dark sides) 3. Specular layer: GGX microsurface reflection 4. Fresnel layer: enhanced environment reflection at grazing angles

Camera Effects Detailed Reference

Prerequisites

• Ray marching fundamentals (ray origin, ray direction)
• Multipass buffers (for accumulation-based DoF)
• Hash functions for stochastic sampling

Thin Lens Model Derivation

A real camera lens focuses light from a focal plane onto the sensor. Points not on the focal plane project to a circle of confusion (CoC) on the sensor.

Circle of Confusion Formula

CoC = |S2 - S1| × A × f / (S1 × (S2 - f))

Where:

• S1 = focal distance (distance to in-focus plane)
• S2 = object distance
• A = aperture diameter
• f = focal length

Simplified for Shaders

CoC ≈ apertureSize × |depth - focalDistance| / depth

Ray-Based Implementation

Instead of computing CoC per pixel, we model the physical process: 1. Choose a random point on the aperture disk → new ray origin 2. The focal point (where the original ray hits the focal plane) stays fixed 3. New ray direction = normalize(focalPoint - newOrigin) 4. Average many such samples → natural bokeh with correct occlusion

Aperture Shape

• Circular: vec2 p = sqrt(r) * vec2(cos(a), sin(a)) — uniform disk
• Polygonal: reject samples outside polygon for hexagonal/octagonal bokeh
• The sqrt(r) is critical for uniform distribution (area-preserving)

Poisson Disk Sampling

Pre-computed 16-point Poisson disk for blur kernels:

const vec2 poissonDisk[16] = vec2[](
    vec2(-0.94201624, -0.39906216), vec2(0.94558609, -0.76890725),
    vec2(-0.09418410, -0.92938870), vec2(0.34495938,  0.29387760),
    vec2(-0.91588581,  0.45771432), vec2(-0.81544232, -0.87912464),
    vec2(-0.38277543,  0.27676845), vec2(0.97484398,  0.75648379),
    vec2(0.44323325, -0.97511554),  vec2(0.53742981, -0.47373420),
    vec2(-0.26496911, -0.41893023), vec2(0.79197514,  0.19090188),
    vec2(-0.24188840,  0.99706507), vec2(-0.81409955,  0.91437590),
    vec2(0.19984126,  0.78641367),  vec2(0.14383161, -0.14100790)
);

Advantages over regular grid: no structured aliasing patterns, better coverage per sample count.

Motion Blur Approaches

Stochastic Time Sampling (Ray Marching)

For each pixel, pick a random time within the shutter interval:

t_sample = iTime + (rand - 0.5) * shutterDuration

Use t_sample for all scene animation. Accumulate multiple frames for convergence.

Velocity Buffer (Post-Process)

1. Render scene + store per-pixel velocity vectors 2. For each pixel, sample along the velocity direction 3. Weight samples by distance from center (triangle filter)

Hybrid

Use temporal accumulation (TAA-style) with per-frame time jitter — converges over frames with no per-frame cost increase.

Film Grain Characteristics

Real film grain properties:

• Luminance-dependent: More visible in shadows, less in highlights
• Temporally varying: Different pattern each frame (use fract(iTime) in hash seed)
• Spatially uncorrelated: Use pixel coordinates in hash, not UV (grain should be screen-resolution)
• Intensity: 0.02-0.05 for subtle, 0.1+ for stylized/vintage look

Cellular Automata & Reaction-Diffusion — Detailed Reference

This document is a detailed supplement to SKILL.md, containing prerequisites, step-by-step explanations, variant details, performance analysis, and complete code examples for combination suggestions.

---

Prerequisites

GLSL Basics

• Uniform variables: iResolution (viewport resolution), iFrame (current frame number), iTime (elapsed time), iMouse (mouse position)
• Texture sampling: texture(iChannel0, uv) samples using UV coordinates (with filtering), texelFetch(iChannel0, ivec2(px), 0) samples at exact integer pixel coordinates
• Multi-buffer feedback architecture: ShaderToy supports Buffer A~D, each buffer can bind itself or other buffers as iChannel input

ShaderToy Multi-Pass Mechanism

Data written by Buffer A → next frame Buffer A reads via iChannel0 self-feedback. This is the core mechanism for inter-frame state persistence. The Image pass handles final visual output.

2D Grid Sampling

• Pixel coordinates fragCoord are floating point, range [0.5, resolution - 0.5]
• UV coordinates = fragCoord / iResolution.xy, range [0, 1]
• texelFetch(iChannel0, ivec2(px), 0) reads the specified pixel exactly (no filtering), suitable for discrete CA
• texture(iChannel0, uv) uses hardware bilinear interpolation, suitable for continuous RD

Basic Vector Math

• normalize(v): normalize a vector
• dot(a, b): dot product
• cross(a, b): cross product
• length(v): vector length

Convolution Kernel Concepts

A 3x3 stencil performs a weighted sum of the center pixel and its 8 neighbors. Different weights produce different effects:

• Laplacian kernel: Detects deviation of the current value from the neighborhood mean (diffusion)
• Gaussian kernel: Blur/smoothing
• Sobel kernel: Edge detection/gradient computation

---

Implementation Steps in Detail

Step 1: Grid State Storage and Self-Feedback

What: Use ShaderToy's Buffer self-read mechanism to persistently store simulation state in a buffer texture. Each frame reads the previous frame's state, computes new state, and writes it back.

Why: GPU shaders are inherently stateless; buffer inter-frame feedback is required for time-step iteration. State is stored in RGBA channels — CA can use a single channel for alive/dead, while RD uses two channels for u and v respectively.

Code:

// Buffer A: read previous frame's own state
// iChannel0 is bound to Buffer A itself (self-feedback)
vec4 prevState = texelFetch(iChannel0, ivec2(fragCoord), 0);

// Can also sample with UV coordinates (supports texture filtering)
vec2 uv = fragCoord / iResolution.xy;
vec4 prevSmooth = texture(iChannel0, uv);

Key points:

• texelFetch performs no filtering, reads a single pixel exactly, suitable for discrete CA
• texture uses hardware bilinear interpolation, blending adjacent pixel values near pixel boundaries, suitable for continuous RD
• The four RGBA channels can store different state variables (e.g., u, v, velocity field components, etc.)

Step 2: Initialization (Noise Seeding)

What: Initialize the grid with pseudo-random noise on the first frame (or first few frames) to provide seeds for the simulation.

Why: Both CA and RD need initial perturbation to start evolution. Different initial conditions produce different final patterns. In practice, seeding is often repeated for the first 2~10 frames, since ShaderToy occasionally skips the first frame.

Code:

// Simple hash noise function
float hash1(float n) {
    return fract(sin(n) * 138.5453123);
}

vec3 hash33(in vec2 p) {
    float n = sin(dot(p, vec2(41, 289)));
    return fract(vec3(2097152, 262144, 32768) * n);
}

// Initialization branch in mainImage
if (iFrame < 2) {
    // CA: random binary initialization
    float f = step(0.9, hash1(fragCoord.x * 13.0 + hash1(fragCoord.y * 71.1)));
    fragColor = vec4(f, 0.0, 0.0, 0.0);
} else if (iFrame < 10) {
    // RD: random continuous value initialization
    vec3 noise = hash33(fragCoord / iResolution.xy + vec2(53, 43) * float(iFrame));
    fragColor = vec4(noise, 1.0);
}

Key points:

• hash1 is a simple pseudo-random number generator based on sin, producing values in [0, 1)
• hash33 generates a 3D random vector from 2D coordinates, used for multi-channel RD initialization
• CA initialization uses step(0.9, ...) to produce approximately 10% density of living cells
• RD initialization uses continuous random values, with iFrame added so each frame seeds differently
• Multi-frame seeding (iFrame < 10) ensures sufficiently rich initial perturbation

Step 3: Neighbor Sampling and Laplacian Computation

What: Perform weighted sampling of the current pixel's 8 (or 4) neighbors, computing the Laplacian or neighbor count.

Why: This is the core of CA/RD — local rules drive state updates through neighbor information. The Laplacian describes how much a point's value deviates from the surrounding average, physically corresponding to diffusion. The nine-point stencil is more accurate and isotropic than a simple cross stencil.

Three Sampling Methods Compared:

Method A Code Details:

// Discrete CA neighbor counting using texelFetch for exact reads
int cell(in ivec2 p) {
    ivec2 r = ivec2(textureSize(iChannel0, 0));
    p = (p + r) % r;  // Wrap-around boundary (toroidal topology), left overflow appears on right
    return (texelFetch(iChannel0, p, 0).x > 0.5) ? 1 : 0;
}

ivec2 px = ivec2(fragCoord);
// Moore neighborhood: sum of 8 neighbors
int k = cell(px + ivec2(-1,-1)) + cell(px + ivec2(0,-1)) + cell(px + ivec2(1,-1))
      + cell(px + ivec2(-1, 0))                          + cell(px + ivec2(1, 0))
      + cell(px + ivec2(-1, 1)) + cell(px + ivec2(0, 1)) + cell(px + ivec2(1, 1));

Method B Code Details:

// Nine-point Laplacian stencil (for RD)
// Weights: diagonal 0.5, cross 1.0, center -6.0 (sum = 0, ensuring Laplacian of a constant field is zero)
vec2 laplacian(vec2 uv) {
    vec2 px = 1.0 / iResolution.xy;
    vec4 P = vec4(px, 0.0, -px.x);
    return
        0.5 * texture(iChannel0, uv - P.xy).xy   // bottom-left
      +       texture(iChannel0, uv - P.zy).xy   // bottom
      + 0.5 * texture(iChannel0, uv - P.wy).xy   // bottom-right
      +       texture(iChannel0, uv - P.xz).xy   // left
      - 6.0 * texture(iChannel0, uv).xy           // center
      +       texture(iChannel0, uv + P.xz).xy   // right
      + 0.5 * texture(iChannel0, uv + P.wy).xy   // top-left
      +       texture(iChannel0, uv + P.zy).xy   // top
      + 0.5 * texture(iChannel0, uv + P.xy).xy;  // top-right
}

Method C Code Details:

// 3x3 weighted blur (Gaussian approximation)
// Weights: diagonal 1, cross 2, center 4, total 16
// Uses vec3 swizzle to cleverly encode 9 offset directions
float blur3x3(vec2 uv) {
    vec3 e = vec3(1, 0, -1);  // e.x=1, e.y=0, e.z=-1
    vec2 px = 1.0 / iResolution.xy;
    float res = 0.0;
    // e.xx=(1,1), e.xz=(1,-1), e.zx=(-1,1), e.zz=(-1,-1) → four diagonals
    res += texture(iChannel0, uv + e.xx * px).x + texture(iChannel0, uv + e.xz * px).x
         + texture(iChannel0, uv + e.zx * px).x + texture(iChannel0, uv + e.zz * px).x;       // ×1
    // e.xy=(1,0), e.yx=(0,1), e.yz=(0,-1), e.zy=(-1,0) → four edges
    res += (texture(iChannel0, uv + e.xy * px).x + texture(iChannel0, uv + e.yx * px).x
          + texture(iChannel0, uv + e.yz * px).x + texture(iChannel0, uv + e.zy * px).x) * 2.; // ×2
    // e.yy=(0,0) → center
    res += texture(iChannel0, uv + e.yy * px).x * 4.;                                          // ×4
    return res / 16.0;
}

Step 4: State Update Rules

What: Apply CA rules or RD differential equations based on neighbor information to compute new state values.

Why: This is the core simulation logic. CA uses discrete decisions (birth/survival/death), RD uses continuous differential equations with Euler integration.

CA Rule Details:

Conway's Game of Life B3/S23 means:

• B3 = Birth when 3 neighbors
• S23 = Survive when 2 or 3 neighbors

int e = cell(px);  // current state (0 or 1)
// Equivalent to: if (k==3) born/survive; else if (k==2 && alive) survive; else die
float f = (((k == 2) && (e == 1)) || (k == 3)) ? 1.0 : 0.0;

Generic Bitmask Rules: Bitmasks can encode arbitrary CA rule sets without modifying logic code. For example:

• B3/S23 → bornset=8 (binary 1000, bit 3), stayset=12 (binary 1100, bits 2,3)
• B36/S23 → bornset=40 (bits 3,5), stayset=12

// stayset/bornset are bitmasks; bit n=1 means triggered when neighbor count is n
float ff = 0.0;
if (currentAlive) {
    ff = ((stayset & (1 << (k - 1))) > 0) ? float(k) : 0.0;  // survive
} else {
    ff = ((bornset & (1 << (k - 1))) > 0) ? 1.0 : 0.0;       // birth
}

RD Gray-Scott Update Details:

Physical meaning of the Gray-Scott equations:

• Du·∇²u: diffusion of u (spatial smoothing)
• -u·v²: reaction consumption (u decreases when u and v meet)
• F·(1-u): replenishment of u (feed, pulling u back toward 1.0)
• Dv·∇²v: diffusion of v
• +u·v²: reaction production (v increases when u and v meet)
• -(F+k)·v: removal of v (combined decay from kill + feed)

float u = prevState.x;
float v = prevState.y;
vec2 Duv = laplacian(uv) * DIFFUSION;  // DIFFUSION = vec2(Du, Dv)
float du = Duv.x - u * v * v + F * (1.0 - u);
float dv = Duv.y + u * v * v - (F + k) * v;
// Forward Euler integration, clamp to prevent numerical instability
fragColor.xy = clamp(vec2(u + du * DT, v + dv * DT), 0.0, 1.0);

Simplified RD Details: This approach doesn't use the standard Gray-Scott equations, but instead uses gradient-driven displacement and random decay to approximate reaction-diffusion behavior. The results are more organic but less controllable.

float avgRD = blur3x3(uv);
vec2 pwr = (1.0 / iResolution.xy) * 1.5;
// Compute gradient (similar to Sobel)
vec2 lap = vec2(
    texture(iChannel0, uv + vec2(pwr.x, 0)).y - texture(iChannel0, uv - vec2(pwr.x, 0)).y,
    texture(iChannel0, uv + vec2(0, pwr.y)).y - texture(iChannel0, uv - vec2(0, pwr.y)).y
);
uv = uv + lap * (1.0 / iResolution.xy) * 3.0;  // Displace sampling point along gradient (diffusion)
float newRD = texture(iChannel0, uv).x + (noise.z - 0.5) * 0.0025 - 0.002;  // Random decay
newRD += dot(texture(iChannel0, uv + (noise.xy - 0.5) / iResolution.xy).xy, vec2(1, -1)) * 0.145;  // Reaction term

Step 5: Visualization and Coloring

What: Map simulation buffer data to visual effects — color mapping, gradient lighting, bump mapping, etc.

Why: Raw simulation data consists of scalar/vector values in 0~1 range, requiring artistic processing to produce appealing visuals. The most common technique is computing the gradient of buffer values to obtain normal information for bump lighting.

Color mapping techniques:

// Basic: nonlinear color separation
// c is a [0,1] value; different pow exponents make RGB channels respond at different rates
float c = 1.0 - texture(iChannel0, uv).y;
vec3 col = pow(vec3(1.5, 1, 1) * c, vec3(1, 4, 12));
// R channel responds linearly, G channel with 4th power (rapid decay in dark areas), B channel with 12th power (blue only at brightest spots)

Gradient normal computation:

// Compute surface normals from scalar field (for bump map lighting)
vec3 normal(vec2 uv) {
    vec3 delta = vec3(1.0 / iResolution.xy, 0.0);
    // Central difference for x and y gradients
    float du = texture(iChannel0, uv + delta.xz).x - texture(iChannel0, uv - delta.xz).x;
    float dv = texture(iChannel0, uv + delta.zy).x - texture(iChannel0, uv - delta.zy).x;
    // z component controls bump intensity (smaller = stronger bumps)
    return normalize(vec3(du, dv, 1.0));
}

Specular highlight effect:

// Produce specular edges via sampling offset
float c2 = 1.0 - texture(iChannel0, uv + 0.5 / iResolution.xy).y;
// c2*c2 - c*c is positive at gradient changes, producing edge highlights
col += vec3(0.36, 0.73, 1.0) * max(c2 * c2 - c * c, 0.0) * 12.0;

Vignette + gamma correction:

// Vignette: darken edges
col *= pow(16.0 * uv.x * uv.y * (1.0 - uv.x) * (1.0 - uv.y), 0.125) * 1.15;
// Fade-in effect
col *= smoothstep(0.0, 1.0, iTime / 2.0);
// Gamma correction (approximately 2.0)
fragColor = vec4(sqrt(min(col, 1.0)), 1.0);

---

Variant Details

Variant 1: Conway's Game of Life (Discrete CA)

Difference from base version: Uses discrete binary state and neighbor counting rules instead of continuous RD equations. This is the most classic cellular automaton, with simple rules that can give rise to extremely complex behavior (gliders, oscillators, still lifes, etc.).

Complete Buffer A code:

int cell(in ivec2 p) {
    ivec2 r = ivec2(textureSize(iChannel0, 0));
    p = (p + r) % r;  // wrap-around boundary
    return (texelFetch(iChannel0, p, 0).x > 0.5) ? 1 : 0;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    ivec2 px = ivec2(fragCoord);

    // Moore neighborhood counting
    int k = cell(px+ivec2(-1,-1)) + cell(px+ivec2(0,-1)) + cell(px+ivec2(1,-1))
          + cell(px+ivec2(-1, 0))                        + cell(px+ivec2(1, 0))
          + cell(px+ivec2(-1, 1)) + cell(px+ivec2(0, 1)) + cell(px+ivec2(1, 1));
    int e = cell(px);

    // B3/S23 rule
    float f = (((k == 2) && (e == 1)) || (k == 3)) ? 1.0 : 0.0;

    // Initialization: approximately 10% random living cells
    if (iFrame < 2) {
        f = step(0.9, fract(sin(fragCoord.x * 13.0 + sin(fragCoord.y * 71.1)) * 138.5));
    }

    fragColor = vec4(f, 0.0, 0.0, 1.0);
}

Adjustment directions:

• Modifying B/S rule numbers can produce completely different behavior
• Increasing initial density (changing the 0.9 in step(0.9, ...)) alters the evolution result
• The .y channel can store "age" for color mapping during visualization

Variant 2: Configurable Rule Set CA (Birth/Survival Bitmask)

Difference from base version: Uses bitmasks to encode arbitrary CA rules, supporting Moore/von Neumann/extended neighborhoods, capable of producing worms, sponges, explosions, and other patterns.

Bitmask encoding explanation:

• BORN_SET = 8 is binary 0b1000, meaning bit 3 is set → B3 (birth when 3 neighbors)
• STAY_SET = 12 is binary 0b1100, meaning bits 2,3 are set → S23 (survive when 2 or 3 neighbors)
• LIVEVAL controls the living cell's state value; when greater than 1, combined with DECIMATE it can produce gradient decay effects
• DECIMATE is the per-frame decay amount, producing a "trailing" effect

Key code:

#define BORN_SET  8        // birth bitmask, 8 = B3 (bit 3 set)
#define STAY_SET  12       // survival bitmask, 12 = S23 (bits 2,3 set)
#define LIVEVAL   2.0      // living cell state value
#define DECIMATE  1.0      // decay value (0=no decay)

// Rule evaluation
float ff = 0.0;
float ev = texelFetch(iChannel0, px, 0).w;
if (ev > 0.5) {
    // Living cell: decay first, then check if survival rule is met
    if (DECIMATE > 0.0) ff = ev - DECIMATE;
    if ((STAY_SET & (1 << (k - 1))) > 0) ff = LIVEVAL;
} else {
    // Dead cell: check if birth rule is met
    ff = ((BORN_SET & (1 << (k - 1))) > 0) ? LIVEVAL : 0.0;
}

Notable rule sets:

• B3/S23 (Conway Life): BORN=8, STAY=12
• B36/S23 (HighLife): BORN=40, STAY=12 — has self-replicators
• B1/S1 (Gnarl): BORN=2, STAY=2 — fractal growth
• B3/S012345678 (Life without death): BORN=8, STAY=511 — only grows, never dies

Variant 3: Separable Gaussian Blur RD (Multi-Buffer Architecture)

Difference from base version: Replaces the single 3x3 Laplacian with separable horizontal/vertical Gaussian blur for the diffusion step, achieving a larger effective diffusion radius with smoother patterns.

Architecture:

• Buffer A: Reaction step (reads Buffer C's blur result as diffusion term)
• Buffer B: Horizontal Gaussian blur (reads Buffer A)
• Buffer C: Vertical Gaussian blur (reads Buffer B)

Why separate:

• A direct NxN kernel requires N² samples
• Separating into horizontal + vertical passes requires N samples each, 2N total
• A 9-tap separable blur = 18 samples ≈ equivalent to an 81-point 9x9 kernel

Buffer B complete code (horizontal blur):

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    vec2 uv = fragCoord / iResolution.xy;
    float h = 1.0 / iResolution.x;
    vec4 sum = vec4(0.0);
    // 9-tap Gaussian weights (approximate normal distribution)
    sum += texture(iChannel0, fract(vec2(uv.x - 4.0*h, uv.y))) * 0.05;
    sum += texture(iChannel0, fract(vec2(uv.x - 3.0*h, uv.y))) * 0.09;
    sum += texture(iChannel0, fract(vec2(uv.x - 2.0*h, uv.y))) * 0.12;
    sum += texture(iChannel0, fract(vec2(uv.x - 1.0*h, uv.y))) * 0.15;
    sum += texture(iChannel0, fract(vec2(uv.x,         uv.y))) * 0.16;
    sum += texture(iChannel0, fract(vec2(uv.x + 1.0*h, uv.y))) * 0.15;
    sum += texture(iChannel0, fract(vec2(uv.x + 2.0*h, uv.y))) * 0.12;
    sum += texture(iChannel0, fract(vec2(uv.x + 3.0*h, uv.y))) * 0.09;
    sum += texture(iChannel0, fract(vec2(uv.x + 4.0*h, uv.y))) * 0.05;
    fragColor = vec4(sum.xyz / 0.98, 1.0);  // 0.98 = weight sum, normalized
}

Buffer C has identical structure but blurs along the y-axis (replace uv.x ± n*h with uv.y ± n*v, where v = 1.0/iResolution.y).

Variant 4: Continuous Differential Operator CA (Vein/Fluid Style)

Difference from base version: Computes curl, divergence, and Laplacian on the grid, combined with multi-step advection loops, producing vein/fluid-like organic patterns that sit between CA and PDE fluid simulation.

Core concepts:

• Curl: Describes the rotational tendency of a field, used to produce vortex effects
• Divergence: Describes the spreading/converging tendency of a field
• Advection: Propagates field values along the velocity field direction

Parameter tuning guide:

• STEPS (10~60): Advection steps; more = smoother but slower
• ts (0.1~0.5): Advection rotation strength, controls vortex intensity
• cs (-3~-1): Curl scaling; negative values produce counter-clockwise rotation
• ls (0.01~0.1): Laplacian scaling, controls diffusion strength
• amp (0.5~2.0): Self-amplification coefficient
• upd (0.2~0.6): Update smoothing coefficient, controls old/new state blend ratio

Key code:

#define STEPS 40
#define ts    0.2
#define cs   -2.0
#define ls    0.05
#define amp   1.0
#define upd   0.4

// Discrete curl and divergence on a 3x3 stencil
// Standard weights: _K0=-20/6 (center), _K1=4/6 (edge), _K2=1/6 (corner)
curl = uv_n.x - uv_s.x - uv_e.y + uv_w.y
     + _D * (uv_nw.x + uv_nw.y + uv_ne.x - uv_ne.y
           + uv_sw.y - uv_sw.x - uv_se.y - uv_se.x);
div  = uv_s.y - uv_n.y - uv_e.x + uv_w.x
     + _D * (uv_nw.x - uv_nw.y - uv_ne.x - uv_ne.y
           + uv_sw.x + uv_sw.y + uv_se.y - uv_se.x);

// Multi-step advection loop
for (int i = 0; i < STEPS; i++) {
    advect(off, vUv, texel, curl, div, lapl, blur);
    offd = rot(offd, ts * curl);  // rotate offset direction
    off += offd;                   // accumulate offset
    ab += blur / float(STEPS);    // accumulate blurred value
}

Variant 5: RD-Driven 3D Surface (Raymarched RD)

Difference from base version: 2D RD results serve as a texture mapped onto a 3D sphere, driving surface displacement and color; the Image pass becomes a full raymarcher.

Implementation points: 1. Buffer A maintains the standard RD simulation unchanged 2. Image pass becomes a raymarching renderer 3. The SDF function maps 3D points to spherical UV, then samples the RD buffer 4. RD values drive surface displacement

Key code:

// Image pass: use RD texture for displacement in the SDF
vec2 map(in vec3 pos) {
    vec3 p = normalize(pos);
    vec2 uv;
    // Spherical parameterization: 3D point → 2D UV
    uv.x = 0.5 + atan(p.z, p.x) / (2.0 * 3.14159);  // longitude [0, 1]
    uv.y = 0.5 - asin(p.y) / 3.14159;                 // latitude [0, 1]

    float y = texture(iChannel0, uv).y;     // read v component from RD buffer
    float displacement = 0.1 * y;            // displacement amount (adjustable scale factor)
    float sd = length(pos) - (2.0 + displacement);  // base sphere SDF + displacement
    return vec2(sd, y);  // return distance and material parameter
}

Extension directions:

• Replace the sphere with a torus, plane, or other base shapes
• Use the two RD channels to separately drive displacement and color
• Add normal perturbation for finer surface detail
• Combine with environment maps for reflection/refraction

---

Performance Optimization In-Depth Analysis

1. texelFetch vs texture

Discrete CA should use texelFetch(iChannel0, ivec2(px), 0) instead of texture():

• Avoids unnecessary texture filtering overhead
• Guarantees pixel-precise reads without floating-point precision causing sampling of adjacent pixels
• For binary states (0/1), any interpolation introduces errors

Continuous RD can use texture() with linear filtering:

• Hardware automatically performs bilinear interpolation
• The interpolation effect is equivalent to additional smoothing/diffusion, which can be advantageous in some cases
• Hardware-accelerated, faster than manual interpolation

2. Separable Blur Instead of Large-Kernel Laplacian

If a large diffusion radius is needed:

• Don't use a larger NxN Laplacian kernel → O(N²) samples
• Do use separable two-pass Gaussian blur (horizontal + vertical) → O(2N) samples
• Implemented through additional buffer passes

Numerical comparison:

3. Multi-Step Sub-Iteration

For RD, you can loop multiple sub-iterations within a single frame using smaller DT, improving convergence speed while maintaining stability:

#define SUBSTEPS 4     // sub-iteration count
#define SUB_DT 0.25    // = DT / SUBSTEPS
for (int i = 0; i < SUBSTEPS; i++) {
    vec2 lap = laplacian9(uv);
    float uvv = u * v * v;
    u += (DU * lap.x - uvv + F * (1.0 - u)) * SUB_DT;
    v += (DV * lap.y + uvv - (F + K) * v) * SUB_DT;
}

Note: In sub-iterations, the Laplacian is only correct when read from the texture on the first step; subsequent steps should recompute the Laplacian based on updated values. However, in practice, the approximation of single-read multi-step integration is often good enough.

4. Reduced-Resolution Simulation

If the target display resolution is high but the pattern's spatial frequency doesn't require 1:1 pixel precision:

• Run the simulation at lower resolution in the buffer (not directly configurable in ShaderToy, but possible in custom engines)
• Use bilinear interpolation upsampling in the Image pass
• Can save 4x~16x computation

5. Avoiding Branches and Conditionals

Use step(), mix(), clamp() instead of if/else for CA rule evaluation to reduce GPU warp divergence:

// Original if/else version:
// if (k==3) f=1.0; else if (k==2 && e==1) f=1.0; else f=0.0;

// Branch-free version:
float f = max(step(abs(float(k) - 3.0), 0.5),
              step(abs(float(k) - 2.0), 0.5) * step(0.5, float(e)));

Explanation:

• step(abs(float(k) - 3.0), 0.5) is 1.0 when k=3, otherwise 0.0
• step(abs(float(k) - 2.0), 0.5) * step(0.5, float(e)) is 1.0 when k=2 and e=1
• max() combines the two conditions

---

Combination Suggestions — Full Details

1. RD + Raymarching (3D Displacement/Shaping)

Map RD results as a heightmap onto 3D surfaces (sphere, plane, torus) and create organic bumpy surfaces through SDF displacement. Suitable for biological organisms, alien terrain, and similar effects.

Complete Image pass example (sphere + RD displacement):

vec2 map(in vec3 pos) {
    vec3 p = normalize(pos);
    vec2 uv;
    uv.x = 0.5 + atan(p.z, p.x) / (2.0 * 3.14159);
    uv.y = 0.5 - asin(p.y) / 3.14159;
    float y = texture(iChannel0, uv).y;
    float displacement = 0.1 * y;
    float sd = length(pos) - (2.0 + displacement);
    return vec2(sd, y);
}

// Use map() in the raymarch loop
// Normals computed via central difference of map()
// Material color based on y value returned by map() for color mapping

2. CA/RD + Particle Systems

Use CA/RD fields as velocity fields or spawn probability fields for particles:

• Particles flow along RD gradients
• New particles spawn at living CA cells
• Produces "living" particle effects

Implementation approach:

• Buffer A: RD/CA simulation
• Buffer B: Particle position storage (each pixel stores one particle's position)
• Image: Visualize particles and/or fields

3. RD + Post-Processing Lighting

In the Image pass, compute normals from RD values → bump mapping → lighting/reflection/refraction. Combined with environment maps (cubemaps), this can produce etched metal surfaces, liquid ripples, and similar effects.

Key techniques:

• Compute gradients from RD scalar field to get normals
• Use Phong/Blinn-Phong lighting model
• Normals used to sample cubemaps for environment reflections
• Multiple color mapping schemes increase visual richness

4. CA + Color Decay Trails

Living cells use high values; after death, values decay each frame (instead of immediately dropping to zero), with different decay rates in RGB channels producing colorful trailing effects. This is the core technique of the Automata X Showcase.

Implementation code example:

// Add decay logic after CA update
vec4 prev = texelFetch(iChannel0, px, 0);
if (f > 0.5) {
    // Living cell: set to high value
    fragColor = vec4(1.0, 1.0, 1.0, 1.0);
} else {
    // Dead cell: different decay rates per channel
    fragColor = vec4(
        prev.x * 0.99,   // R decays slowly → longest red trail
        prev.y * 0.95,   // G decays moderately
        prev.z * 0.90,   // B decays fast → shortest blue trail
        1.0
    );
}

5. RD + Domain Warping

Apply vortex warp or spiral zoom domain transforms to the RD sampling UV before computing, causing the diffusion field itself to be distorted, producing spiral and vortex-like organic patterns. Flexi's Expansive RD uses this technique.

Implementation code example:

// Apply domain transform to UV before RD update
vec2 warpedUV = uv;
// Vortex warp
float angle = length(uv - 0.5) * 3.14159 * 2.0;
float s = sin(angle * 0.1);
float c = cos(angle * 0.1);
warpedUV = (warpedUV - 0.5) * mat2(c, -s, s, c) + 0.5;

// Sample state using transformed UV
vec2 state = texture(iChannel0, warpedUV).xy;
// Then proceed with normal RD computation...

SDF Tricks Detailed Reference

Prerequisites

• Understanding of signed distance fields and ray marching
• Basic SDF primitives and boolean operations
• FBM / procedural noise fundamentals

Lipschitz Condition and FBM Detail

An SDF must satisfy the Lipschitz condition: |f(a) - f(b)| ≤ |a - b| (gradient magnitude ≤ 1). This guarantees that stepping by the SDF value is always safe — no surface exists within that radius.

When adding FBM noise to an SDF, the noise derivatives can violate Lipschitz:

• Raw noise amplitude of 0.1 with frequency 20 has gradient ~2.0, breaking the condition
• This causes ray marching to overshoot, creating holes and artifacts

Solutions: 1. Amplitude limiting: Keep amplitude × frequency < 1.0 across all octaves 2. Distance fade: d += amp * fbm(p * freq) * smoothstep(fadeStart, 0.0, d) — detail only appears near the surface where overshoot distance is small 3. Step size reduction: Multiply ray step by 0.5-0.7, trading speed for stability

Bounding Volume Strategies

Hierarchical Bounding

For scenes with N objects, test bounding volumes in order of increasing cost:

Level 1: Scene bounding sphere (1 evaluation)
Level 2: Object group bounds (few evaluations)
Level 3: Individual object SDF (full cost)

Spatial Partitioning

For repeating structures, combine domain repetition with bounds:

float map(vec3 p) {
    vec3 q = mod(p + 2.0, 4.0) - 2.0;  // repeat every 4 units
    // Only evaluate detail if within local bounding sphere
    float bound = length(q) - 1.5;
    if (bound > 0.2) return bound;
    return detailedSDF(q);
}

Binary Search Convergence

After N iterations of binary search, the position error is initialStep / 2^N:

• 4 iterations: 1/16 of initial step size
• 6 iterations: 1/64 of initial step size (sub-pixel at typical resolutions)
• 8 iterations: 1/256 (overkill for most uses)

6 iterations is the practical sweet spot — gives sub-pixel precision without wasting GPU cycles.

XOR Operation Mathematics

opXor(a, b) = max(min(a, b), -max(a, b))

This is equivalent to: union(a, b) AND NOT intersection(a, b) — the symmetric difference. Geometry exists where exactly one shape is present but not both. Useful for creating lattice structures and interlocking patterns.

Interior SDF Pattern Techniques

When the camera is inside an SDF (d < 0), the negative distance still gives useful information:

• abs(d) gives distance to nearest surface from inside
• Combine with repeating patterns using fract() to create infinite interior structures
• Use max(outerSDF, innerSDF) to confine interior patterns within the outer shell

Advanced Texture Mapping Detailed Reference

Prerequisites

• Screen-space derivatives (dFdx, dFdy)
• textureGrad() function usage
• Basic ray marching

Triplanar vs Biplanar Cost Analysis

Modern GPUs are typically bandwidth-limited rather than ALU-limited, making biplanar the better default choice.

Weight Remapping Mathematics

The biplanar weight formula clamp((w - 0.5773) / (1.0 - 0.5773), 0, 1) ensures:

• At normals aligned with one axis: weight = 1.0 (clean projection)
• At 45° diagonals where 2 axes are equal: smooth transition
• At the cube diagonal (1/√3 ≈ 0.5773): weight = 0.0, but this is the point where the third (discarded) projection would be needed — biplanar's approximation error is maximal here but visually acceptable

Gradient Propagation

Using textureGrad() instead of texture() is essential because: 1. Axis selection (ma, me) creates UV discontinuities at projection boundaries 2. Hardware texture() computes mip from implicit derivatives, which spike at discontinuities → visible seams 3. textureGrad() with manually propagated dFdx(p), dFdy(p) bypasses this, keeping gradients smooth across boundaries

Ray Differential Mathematics

Problem Statement

In rasterization, dFdx/dFdy of texture coordinates work naturally because adjacent pixels map to nearby surface points. In ray marching, adjacent pixels may hit completely different objects → broken mip selection.

Solution: Tangent Plane Intersection

Given:

• Primary ray hits surface at pos with normal nor
• Neighbor pixel ray rd_neighbor originates from ro_neighbor

The neighbor ray's intersection with the tangent plane at pos:

t_neighbor = dot(pos - ro_neighbor, nor) / dot(rd_neighbor, nor)
pos_neighbor = ro_neighbor + rd_neighbor * t_neighbor

The difference pos_neighbor - pos gives the world-space footprint of one pixel at the hit point.

For Perspective Cameras (Common Case)

ro is the same for all pixels, only rd varies:
dposdx = t * (rdx * dot(rd, nor) / dot(rdx, nor) - rd)
dposdy = t * (rdy * dot(rd, nor) / dot(rdy, nor) - rd)

Where rdx = rd + dFdx(rd) and rdy = rd + dFdy(rd).

Chain Rule for Texture Coordinates

If texture mapping function is uv = f(pos):

duvdx = Jacobian(f) × dposdx
duvdy = Jacobian(f) × dposdy

For simple planar mapping uv = pos.xz:

duvdx = dposdx.xz
duvdy = dposdy.xz

Texture Repetition Theory

Why Tiling is Visible

Human vision excels at detecting: 1. Periodic patterns: Regular grid alignment 2. Unique features: Distinctive spots/marks that repeat identically 3. Phase alignment: All tiles start at the same phase

Breaking Repetition

Each method targets different cues:

• Random offset (Method A): Breaks phase alignment, 4 fetches
• Voronoi blend: Breaks grid structure entirely, 9 fetches (expensive)
• Virtual pattern (Method B): Breaks unique features cheaply, 2 fetches

Method B is preferred for real-time use — the low-frequency index variation is cache-friendly and the two texture fetches share locality.

WebGL2 Pitfalls Reference

This is a reference document for the webgl-pitfalls technique.

Complete Error Message Reference

Type Mismatch Examples

// ERROR: terrainM expects vec2, passing vec3
float calcAO(vec3 pos, vec3 nor) {
    float d = terrainM(pos + h * nor);  // Wrong: pos + h*nor is vec3
}
// FIX: Extract xz components
float calcAO(vec3 pos, vec3 nor) {
    float d = terrainM(pos.xz + h * nor.xz);  // Correct: vec2
}

// ERROR: can't access .z on vec2
vec2 uv = vec2(1.0, 2.0);
float z = uv.z;  // Wrong: vec2 has no .z
// FIX: use proper swizzle or conversion
float z = uv.y;  // Or if you need third component, use vec3

GLSL ES 3.0 Specific Notes

• All declared uniform variables must be used in shader code, otherwise compiler may optimize them away
• When gl.getUniformLocation() returns null, setting that uniform triggers INVALID_OPERATION
• Loop counters must be deterministic at runtime — avoid compile-time constant folding issues

Anti-Aliasing Techniques

Use Cases

• Eliminating jagged edges (staircase artifacts) in ray-marched or SDF-rendered scenes
• Smooth 2D SDF shape rendering
• Post-process edge smoothing for any shader output
• Temporal smoothing for noise reduction

Core Principles

Anti-aliasing in shaders differs from rasterization pipelines. Without hardware MSAA on procedural geometry, we rely on analytical or post-process approaches.

Techniques

1. Supersampling (SSAA) for Ray Marching

Render multiple sub-pixel samples and average:

#define AA 2  // 1=off, 2=4x, 3=9x
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    vec3 totalColor = vec3(0.0);
    for (int m = 0; m < AA; m++)
    for (int n = 0; n < AA; n++) {
        vec2 offset = vec2(float(m), float(n)) / float(AA) - 0.5;
        vec2 uv = (2.0 * (fragCoord + offset) - iResolution.xy) / iResolution.y;
        vec3 col = render(uv);
        totalColor += col;
    }
    fragColor = vec4(totalColor / float(AA * AA), 1.0);
}

Cost: AA^2 × full render. Use AA=2 for quality, AA=1 for development.

2. SDF Analytical Anti-Aliasing

For 2D SDF shapes, use pixel width to compute smooth edges:

float d = sdShape(uv);
float fw = fwidth(d);  // screen-space derivative of SDF
float alpha = smoothstep(fw, -fw, d);  // smooth edge over exactly 1 pixel

// Alternative: manual pixel width for more control
float pixelWidth = 2.0 / iResolution.y;  // approximate pixel size in UV space
float alpha2 = smoothstep(pixelWidth, -pixelWidth, d);

For 3D SDF scenes, apply anti-aliasing at the edge of geometry:

// After ray marching, at the surface:
float edgeFade = 1.0 - smoothstep(0.0, 0.01 * t, lastSdfValue);
// t = ray distance — scales threshold with distance for consistent edge width

3. Temporal Anti-Aliasing (TAA) Basics

Blend current frame with previous frame using a multipass buffer:

// Buffer A: render with sub-pixel jitter
vec2 jitter = (hash22(vec2(iFrame)) - 0.5) / iResolution.xy;
vec2 uv = (fragCoord + jitter) / iResolution.xy;
vec3 currentColor = render(uv);

// Buffer A output: store current render
fragColor = vec4(currentColor, 1.0);

// Image shader: blend with history
vec3 current = texture(iChannel0, fragCoord / iResolution.xy).rgb;  // this frame
vec3 history = texture(iChannel1, fragCoord / iResolution.xy).rgb;  // previous frame
float blend = 0.9;  // higher = smoother but more ghosting
fragColor = vec4(mix(current, history, blend), 1.0);

Note: Full TAA also needs motion vectors and neighborhood clamping to avoid ghosting.

4. FXAA (Fast Approximate Anti-Aliasing)

Simplified post-process edge detection and smoothing:

vec3 fxaa(sampler2D tex, vec2 uv, vec2 texelSize) {
    // Sample center and 4 neighbors
    vec3 rgbM = texture(tex, uv).rgb;
    vec3 rgbN = texture(tex, uv + vec2(0.0, texelSize.y)).rgb;
    vec3 rgbS = texture(tex, uv - vec2(0.0, texelSize.y)).rgb;
    vec3 rgbE = texture(tex, uv + vec2(texelSize.x, 0.0)).rgb;
    vec3 rgbW = texture(tex, uv - vec2(texelSize.x, 0.0)).rgb;

    // Luma for edge detection
    vec3 lumaCoeff = vec3(0.299, 0.587, 0.114);
    float lumaN = dot(rgbN, lumaCoeff);
    float lumaS = dot(rgbS, lumaCoeff);
    float lumaE = dot(rgbE, lumaCoeff);
    float lumaW = dot(rgbW, lumaCoeff);
    float lumaM = dot(rgbM, lumaCoeff);

    float lumaMin = min(lumaM, min(min(lumaN, lumaS), min(lumaE, lumaW)));
    float lumaMax = max(lumaM, max(max(lumaN, lumaS), max(lumaE, lumaW)));
    float lumaRange = lumaMax - lumaMin;

    // Skip if edge contrast is low
    if (lumaRange < max(0.0312, lumaMax * 0.125)) return rgbM;

    // Blend along edge direction
    vec2 dir;
    dir.x = -((lumaN + lumaS) - 2.0 * lumaM);
    dir.y = ((lumaE + lumaW) - 2.0 * lumaM);
    float dirReduce = max(lumaRange * 0.25, 1.0 / 128.0);
    float rcpDirMin = 1.0 / (min(abs(dir.x), abs(dir.y)) + dirReduce);
    dir = clamp(dir * rcpDirMin, -8.0, 8.0) * texelSize;

    vec3 rgbA = 0.5 * (texture(tex, uv + dir * (1.0/3.0 - 0.5)).rgb +
                        texture(tex, uv + dir * (2.0/3.0 - 0.5)).rgb);
    return rgbA;
}

Choosing the Right Approach

→ For deeper details, see reference/anti-aliasing.md